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1.0  SUMMARY 


Due  to  the  rising  cost  of  launching  payloads  into  outer  space  and  to  address  the  problems 
of  aging  spacecraft,  engineers  are  designing  new  vehicles  that  will  reduce  turnaround  time  and 
cost  per  launch.  This  activity  has  led  to  a  large  increase  in  research  effors  conducted  in  the  field 
of  reusable  launch  vehicles  (RLV).  Vehicle  design  is  an  inherently  non-linear,  multi  physics 
based  problem  involving  continuous,  mixed,  and  integer  optimization  variables.  With  the  wide 
presence  of  randomness  in  these  variables,  including  epistemic,  aleatory,  and  model-form 
uncertainties,  uncertainty  in  the  problem  must  be  accounted  for  in  order  to  accurately  quantify 
the  systems  probability  of  success. 

This  research  provides  valuable  tools  to  the  RLV  community  by  incorporating  risk- 
minimization  into  the  design  of  a  vehicle.  This  has  been  accomplished  by  incorporating 
uncertainty  quantification  related  to  the  aeroelastic  and  structural  integrity  of  a  launch  vehicle, 
with  a  focus  on  flutter  uncertainties. 

In  this  research  a  trajectory  optimization  was  analyzed,  from  which  critical  flight  points 
along  the  mission  path  were  identified.  To  further  investigate  these  critical  points,  a  finite 
element  structural  and  aeroelastic  model  was  created  to  represent  the  RLV  wing  configuration. 

To  ensure  an  accurate  wing  model,  the  model  was  validated  based  on  a  frequency  and  mode 
shape  comparison  analysis. 


1 


Aleatory  uncertainty,  epistemic  uncertainty,  and  model-form  uncertainty  quantification 
techniques  were  explored  for  their  suitability.  To  minimize  risk  associated  with  the  variabilities 
in  aeroelastic  flutter,  a  reliability-based  design  optimization  is  considered  in  this  investigation. 

Three  analyses  were  conducted  on  the  FAST  configuration  F  wing.  The  first  analysis 
includes  aleatory  uncertainty  in  structural  geometry.  The  second  investigates  epistemic 
uncertainty  quantification  incorporating  evidence  theory,  where  uncertainties  in  atmospheric 
conditions  and  composite  materials  were  quantified.  Finally,  a  gradient-based  reliability  design 
optimization  was  investigated,  resulting  in  a  weight  reduction  and  increased  structural  reliability 
of  the  RLV’s  wing. 

The  technical  contents  of  this  report  are  organized  in  the  following  manner.  Chapter  2 
presents  the  background  information  for  an  RLV.  Chapter  3  describes  the  mission  trajectory 
analysis  that  served  as  the  basis  for  defining  critical  aerodynamic  loads.  Chapter  4  discusses  the 
development  of  the  finite  element  model  for  the  tip-tail  FAST  configuration  F  wing.  Chapter  5 
demonstrates  techniques  to  quantify  aleatory,  epistemic  and  model-form  uncertainty.  Chapter  6 
investigates  reliability-based  design  optimization.  Chapter  7  includes  results  for  aleatory  and 
epistemic  uncertainty  analyses  of  the  FAST  vehicle’s  wing  as  well  as  a  reliability-based  design 
optimization.  The  report  then  closes  with  Chapter  8  with  concluding  remarks. 
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2.0  INTRODUCTION 


The  United  States  now  has  an  operating  mixed  fleet  of  space  vehicles  comprised  of 
reusable  space  shuttle  orbiters  and  expendable  launch  vehicles  (ELVs).  To  reduce  operational 
costs,  launch  vehicles  must  be  designed  to  be  reusable.  In  recent  years  the  United  States  Air 
Force  (USAF)  and  NASA  have  conducted  ongoing  research  in  the  design  of  reusable  launch 
vehicles  (RFV)  [1  -  6].  The  US  Air  Force  has  been  developing  a  reusable  launch  vehicle 
variously  called  Micro-X,  Hot  Eagle,  and  now  FAST  (Future  Responsive  Access  to  Space 
Technology).  Configuration  design  activities  for  the  FAST  launch  vehicle  demonstrator  have 
focused  on  configurations  that  have  vertical  tails  on  the  wing  tips  as  seen  in  Figure  1.  This  project 
seeks  to  develop  and  mature  technologies  that  will  allow  development  of  future  military  launch 
vehicles  that  can  land  to  launch  another  day. 


Figure  1:  FAST  Vehicle  Tip-Tail  Design 
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Future  space  access  vehicle  designs  employ  composite  materials,  thennal  protection 
systems,  health  monitoring  hardware,  and  propellant  tanks  to  make  the  vehicles  fully  reusable 
and  reliable.  With  many  integrated  subsystems,  development  of  future  vehicles  with  high 
reliability  and  confidence  requires  investigation  of  the  interaction  of  multi  physics  behavior  as 
well  as  optimization  and  quantification  of  the  risk  associated  with  each  subsystem. 

This  technical  effort  with  AFRL  addresses  two  aspects  of  future  vehicle  design:  First  the 
work  incorporates  the  vehicle  trajectory  by  meeting  the  payload  requirements  at  the  systems 
level.  Second  at  the  subsystem  level,  this  assessment  conducts  an  aeroelastic  flutter  analysis 
which  is  used  as  a  baseline  to  quantify  and  mitigate  the  risk  in  the  aeroelastic  flutter  dynamic 
pressure.  Flutter  is  an  important  characteristic  for  uncertainty  quantification  since  this  type  of 
configuration  is  more  apt  to  flutter  at  lower  velocities  which  occur  during  the  vehicles  trajectory. 
The  important  goal  of  this  assessment  is  to  develop  a  toolbox  which  can  be  used  to  easily 
quantify  different  types  of  uncertainties.  A  spectrum  of  uncertainty  quantification  procedures 
must  be  incorporated  in  the  design  toolbox  to  consider  probability  distributions,  parameter 
intervals,  expert  opinions,  spatial  and  time-dependent  variances,  and  limited  experimental  data  in 
defining  the  subcomponent  models.  Through  the  use  of  advanced  stochastic  expansion 
techniques,  the  designer  is  supplied  with  the  risks  associated  with,  each  of  the  selected 
alternatives.  Some  of  the  metrics  of  measure  are  the  probability  of  failure,  reliability,  confidence 
bounds,  belief  plausibility,  and  sensitivity  factors. 
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3.0  MISSION  GOALS 


The  RLV  mission  goals  are  to  examine  how  uncertainty  propagates  through  system  level 
traits  to  the  subsystem  level.  The  first  step  of  this  process  is  to  identify  a  trajectory  that  the 
vehicle  will  follow,  to  complete  a  desired  mission.  In  this  chapter  a  trajectory  optimization  code 
was  selected  and  executed  based  on  the  mission  requirements  of  a  reusable  launch  vehicle. 
Critical  points  along  the  trajectory  were  identified  and  transmitted  into  more  detailed  FEA 
simulations  in  the  following  chapters. 

3.1  Trajectory  Optimization 

Trajectory  is  the  path  an  object  follows  in  order  to  reach  a  final  destination.  There  are 
many  different  types  of  trajectories  the  FAST  vehicle  can  travel  [7,  8].  The  trajectories  can  range 
from  a  downrange  mission,  where  the  vehicle  launches  from  one  site  and  lands  at  another,  to  a 
rocket-back  and  return  to  launch  site.  To  begin  vehicle  design,  a  trajectory  analysis  is  needed  to 
determine  what  kind  of  flight  loads  the  RLV  will  encounter  during  the  mission.  There  are  many 
codes  available  to  solve  trajectory  optimization.  Several  codes  include  Optimal  Trajectories  by 
Implicit  Simulation  (OTIS),  Program  to  Optimize  Simulated  Trajectories  (POST),  Graphical 
Environment  for  Simulation  and  Optimization/Aerospace  Trajectory  Optimization  Software 
(GESOP/ASTOS),  and  General  Trajectory  Simulation  (GTS).  This  research  will  focus  on  the 
most  widely  used  and  accepted  trajectory  optimization  codes  POST  and  OTIS. 
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3.1.1  POST  and  OTIS  Comparison 


Two  trajectory  optimization  codes  commonly  accepted  in  the  aerospace  community  are 
POST  and  OTIS.  POST  and  OTIS  represent  the  structural  models  of  the  launch  vehicle  with  a 
generalized  point  mass.  The  programs  are  able  to  target  and  optimize  point  mass  trajectories  for 
an  assortment  of  powered  or  unpowered  vehicles.  In  a  recent  journal  paper  these  two  trajectory 
optimization  codes  were  compared  [9].  Physical  and  natural  phenomenons  are  modeled  and 
calculated  the  same  in  both  of  the  optimization  programs.  In  theory  each  program  should 
produce  the  same  results  given  that  they  are  both  searching  for  an  optimum.  Any  small 
discrepancies  in  the  two  optimization  codes  outputs  can  be  attributed  by  the  different 
mathematical  optimizers  and  algorithms  used  during  the  analyses. 

POST  employs  a  gradient-based  approach  by  integrating  each  point  in  the  trajectory 
equation.  During  the  analysis  each  point  in  the  trajectory  must  satisfy  all  physical  laws  to  move 
to  the  next  point  in  the  trajectory.  OTIS  has  two  methods  to  solve  trajectory  optimization,  like 
POST  it  also  contains  a  gradient-based  approach  as  well  as  a  collocation  method.  OTIS’s 
collocation  method  allows  more  flexibility  by  giving  the  user  the  ability  to  choose  the  trajectory 
nodes  and  node  spacing  in  the  program.  In  theory  an  infinite  amount  of  degrees  of  freedom  could 
be  examined,  rather  than  the  program  choosing  the  points  for  the  user.  One  drawback  to  this 
technique  may  not  always  produce  an  actual  physical  solution  to  the  problem. 

Both  programs  have  a  user  interface  where  the  optimization  parameters  are  input.  OTIS 
requires  the  user  to  define  the  continuity,  that  can  lead  to  problems  with  inexperienced  users. 
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POST  automatically  determines  the  continuity  making  it  a  more  applicable  program  to  optimize 
the  FAST  vehicle’s  trajectory. 

Although  the  two  programs  seem  very  similar,  POST  appears  to  be  more  user  friendly 
and  straightforward.  POST  II  was  developed  from  the  original  POST  software  starting  in  1995 
and  was  written  in  FORT  AN  77  and  C  languages  [10]. 

POST  II’s  optimization  procedure  consists  of: 

Inputting  the  data  for  all  the  phases 
Initialize  equations  of  motion 

Propagating  the  trajectory  until  interrupted  by  the  occurrence  of  the  user-specified 
condition  for  the  next  event 
Restart  equations  of  motion  with  the  new  inputs 
Repeat  until  final  event  and  the  post  process  the  data 


POST  II  can  perform  optimization  targeting  with  or  without  inequality  constraints, 
unconstraint  optimization,  and  constrained  (equal  and/or  inequality  constraints)  optimization. 

The  program  employs  finite  differencing  to  calculate  the  sensitivities  of  the  optimization  variable 
with  respect  to  the  constraints.  The  user  selects  the  integration  method  and  initial  step  size.  The 
selection  of  the  integration  and  step  size  are  an  important  optimization  consideration  because  it 
can  make  drastic  impact  on  the  computation  time.  POST  II  provides  an  automatic  perturbation 
step  size  controller,  although  an  initial  guess  is  needed,  it  should  make  the  loss  of  accuracy  much 
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less  than  if  a  wrong  perturbation  variable  was  chosen,  which  can  lead  to  introducing  more 
uncertainty  in  the  problem. 

3.2  FAST  Trajectory  Analysis 

In  this  research  a  trajectory  optimization  using  POST  was  conducted  on  a  nominal 
rocket-back  mission.  This  type  of  mission  requires  the  RLV  to  launch  in  a  vertical  position. 
Upon  completing  the  mission  the  vehicle  returns  to  the  launch  point  and  lands  horizontally  like  a 
conventional  aircraft,  ready  to  be  refueled  and  launched  again  in  a  short  turnaround. 

Figure  2  depicts  the  trajectory  plot  generated  by  POST  II  of  the  RLV’s  trajectory.  This 
trajectory  will  supply  essential  information  to  complete  finite  element  analysis  at  the  critical 
points.  A  user  specifies  how  the  vehicle  performs  and  the  boundary  conditions  for  a  desired 
mission.  POST  II  then  solves  the  optimal  trajectory  based  on  the  constraints  and  vehicle 
perfonnance  data.  Boundary  conditions  include  Mach  number,  final  velocity,  altitudes,  and 
launch  angles. 


8 


,10* 


Figure  2:  RLV  Rocket  Back  Trajectory  Analysis 


The  trajectory  the  RLV  travels  is  explained.  The  vehicle  starts  from  a  vertical  launch  pad  and 
launches.  For  a  short  duration,  the  vehicle  continues  vertically  to  reach  the  desired  altitude  and  Mach 
number  that  defines  the  staging  point.  During  this  time  the  vehicle  experiences  heavy  aerodynamic 
loading.  For  an  operational,  multi-stage  launch  vehicle  additional  stages  would  be  located  on  the 
reusable  booster.  This  would  allow  the  vehicles  to  separate  and  continue  on  their  desired  mission. 

The  propellant  capacities  for  ascent  propellant  and  rocket-back  propellant  were  allocated  to  allow  the 
vehicle  to  return  to  the  launch  site.  After  the  launch  vehicle  and  payload  separates,  the  vehicles  coast 
for  a  short  time  allowing  them  to  have  enough  clearance  to  cause  no  damage.  At  this  time  the 
reusable  booster  will  then  pitch-over  to  an  optimum  altitude  initiate  for  the  reentry  portion  of  the 
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trajectory.  At  the  apex  of  the  trajectory  the  vehicle  changes  orientation  to  allow  for  an  optimal  angle 
of  attack  (AO A)  to  reduce  the  temperatures  the  vehicle  will  encounter  while  entering  the  Earth’s 
atmosphere.  The  vehicle  then  transitions  to  a  maximum  range  glide  altitude  during  the  landing 
stage  as  it  approaches  the  runway  of  the  landing  site. 

3.3  FAST  Critical  Mission  Points 

Since  the  trajectory  optimization  program  only  explores  the  design  space  with  a  point 
mass,  the  perfonnance  analysis  information  obtained  cannot  be  assumed  to  be  exact.  To  quantify 
uncertainty  of  the  vehicles  performance,  a  higher  fidelity  analysis  must  be  utilized  using  finite 
element  analysis  (FEA)  approaches.  Executing  finite  element  analysis  on  the  entire  trajectory 
would  be  extremely  computationally  expensive.  This  calls  for  critical  points  along  the  trajectory 
to  be  recognized  and  analyzed.  Critical  points  along  the  trajectory  were  identified  during  the 
trajectory,  including:  lift  off,  separation  of  the  booster  from  the  payload,  the  rocket  back  phase, 
reentry  phase,  as  well  as  landing.  At  each  of  these  critical  points  an  assortment  of  finite  element 
analysis  techniques  can  be  implemented  depending  on  what  scenario  the  vehicle  is  encountering, 
resulting  in  more  accurate  information  for  that  specific  flight  point.  For  example  a  flutter 
analysis  may  be  necessary  at  reentry  or  launching  due  to  the  vehicles  experience  in  high  dynamic 
pressures,  knowing  flutter  can  be  affected  by  dynamic  pressure.  In  contrary  a  maneuverability 
analysis  may  be  conducted  during  the  booster  separation  point  or  apex  where  the  vehicle  is 
making  precise  maneuvers,  which  may  not  be  physically  possible  although  the  optimization 
program  provides  these  as  potential  paths.  These  critical  points  will  be  analyzed  using 
NASTRAN  and  ASTROS,  finite  element  modeling  packages  that  include  aeroelastic  analysis. 
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4.0  REUSABLE  LAUNCH  VEHICLE  FINITE  ELEMENT 

MODEL  DEVELOPMENT 

To  fully  analyze  the  reusable  launch  vehicle,  representative  models  and  simulations  must 
be  developed  and  implemented  as  tools  in  the  design  process.  This  chapter  demonstrates  the 
modeling  techniques  used  to  capture  the  RLV’s  wing  structure  as  well  as  incorporating 
aeroelastic  flutter  analysis. 

The  RLV  will  have  a  tip-tail  configuration  presented  in  Figure  1.  These  configurations 
have  advantages  related  to  aerodynamic  effectiveness  of  the  tails  during  high  AOA  reentry,  since 
they  are  not  shadowed  by  the  body.  They  also  hold  operability  benefits  in  that  the  tails  are  away  from 
the  aft  end  of  the  vehicle  and  thus  allow  easier  access  to  the  engines.  The  tip-tail  wings  also  offer  the 
opportunity  to  mount  a  payload  system  on  top  of  the  vehicle  without  interference  from  the  vertical 
tails. 

4.1  Structural  Model  Development 

It  is  important  to  be  able  to  identify  and  quantify  modeling  errors  introduced  when 
different  fidelity  simulations  are  utilized.  While  high  fidelity  analyses  generally  produce 
solutions  most  reflective  upon  the  actual  scenario,  the  computational  cost  associated  with  these 
simulations  can  become  restrictive  within  a  design  environment.  Even  for  relatively  simple 
designs,  a  large  number  of  function  evaluations,  which  in  a  simulation-based  environment 
directly  correlates  to  the  number  of  simulations  required  and  necessary  for  risk  quantification 
and  design  optimization,  are  needed.  Thus,  it  is  important  to  be  able  to  utilize  lower  fidelity 
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simulations  where  less  accurate  information  is  required,  and  only  utilize  high  fidelity  simulations 
in  aspects  of  the  design  where  more  exact  information  is  needed,  such  as  using  lower  fidelity 
simulations  in  the  initial  exploration  of  the  design  space,  and  then  utilizing  the  higher  fidelity 
models  as  the  design  space  is  reduced. 

4.1.1  Reusable  Launch  Vehicle  Wing  FEA  Structural  Model 

Figure  3  is  a  high  fidelity  structural  model  developed  by  a  highly  regarded  aircraft  design 
company.  This  model  was  used  as  a  baseline  to  create  the  wing  structural  model  that  will  be  used 
in  the  uncertainty  quantification  design  process. 


k 


Figure  3:  High  Fidelity  RLV  Structural  Model 


The  RLV’s  wing  model  was  modeled  based  on  the  parameters  found  in  Table  1.  In  this 
model,  membrane  elements  with  bending  capability  are  used  to  represent  the  wing  skins,  shear 
panels  represent  the  spars  and  ribs,  and  rod  elements  represent  the  posts.  The  center  of  gravity  is 
maintained  by  adding  non  structured  mass  along  the  line  of  center  of  gravity  axis. 
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Table  1:  RLV  Wing  Parameters 


Span 

21.8  ft 

Root  Chord 

18.2  ft 

Tip  Chord 

6.1  ft 

Mean  Chord 

13.1  ft 

Aspect  Ratio 

1.8 

Taper  Ratio 

0.333 

Area  Sref 

265  ft2 

Sweep 

48  deg 

Dry  Weight 

800  lb 

The  structural  model  was  developed  in  NASTRAN  using  the  following  elements. 

Skins:  QUAD-4/TRIA-3  Elements 

Spars:  Shear  Elements 

Ribs:  Shear  Elements 

Spar  Caps:  Bar  Elements 

Figure  4  is  an  image  of  the  RLV’s  wing  obtained  from  full  RLV  structural  model  which 
will  be  referred  to  as  the  original  model.  The  elements  that  construct  the  model  in  NASTAN  can 
be  found  below  Figure  4. 
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Figure  4:  Original  Configuration  F  Model 


NUMBER 

OF 

GRID 

POINTS 

465 

NUMBER 

OF 

CONM2 

ELEMENTS  = 

12 

NUMBER 

OF 

CQUAD4 

ELEMENTS  = 

204 

NUMBER 

OF 

CROD 

ELEMENTS  = 

605 

NUMBER 

OF 

CSHEAR 

ELEMENTS  = 

174 

NUMBER 

OF 

CTRIA3 

ELEMENTS  = 

28 

NUMBER 

OF 

RBAR 

ELEMENTS  = 

133 

NUMBER 

OF 

RBE1 

ELEMENTS  = 

5 

To  reduce  the  complexity  of  the  original  model,  elements  were  removed  to  produce  a 


more  streamline  model  seen  in  Figure  5.  The  modified  model  of  the  RLV’s  wing  will  reduce  the 


computational  time  when  employing  uncertainty  quantification  methods.  Once  again  the 


elements  used  in  the  modified  model  can  be  found  below  the  Figure  5. 
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Figure  5:  Modified  Configuration  F  Model 
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OF 

GRID 

POINTS 

91 
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OF 

CONM2 
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5 

NUMBER 
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CQUAD4 
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54 
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OF 

CROD 

ELEMENTS  = 

201 
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OF 

CSHEAR 
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72 

NUMBER 

OF 

CTRIA3 

ELEMENTS  = 

6 

This  modified  model  is  capable  of  static  and  dynamic  aeroelasticty  simulations  with 
possible  analysis  including: 

-  Static  Analysis: 

-  Divergent  Dynamic  Pressure 

-  Trim  Angle  of  Attack 

-  Control  Surface  Reversal  Pressure 

-  Structural  Displacement 

-  Dynamic  Analysis: 

-  Flutter  Velocity 

-  Low-frequency  modes  and  mode  shape 
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4.2  Structural  Model  Validation 


When  a  model  is  modified  from  its  original  configuration  it  is  necessary  to  know  how  the 
two  models  compare.  To  be  able  to  use  the  modified  model  with  confidence  a  validation  study 
must  be  conducted.  Two  methods  of  validation  were  demonstrated  in  the  following  sections.  The 
first  validation  implemented  was  a  natural  frequency  comparison,  where  the  two  wing  models 
first  five  natural  frequencies  were  analyzed.  The  second  validation  uses  the  wing  models  mode 
shapes  to  investigate  a  modal  assurance  criterion  analysis.  Although  both  types  of  validation  can 
be  used  separately  a  more  conclusive  comparison  can  be  made  when  the  two  methods  are 
coupled. 


4.2.1  Frequency  Comparison  Analysis 


The  first  validation  analysis  compared  the  natural  frequencies  between  the  two  models 
described  in  section  4.1.1.  The  natural  frequencies  were  established  using  NASTRAN’s  dynamic 
analysis.  Both  models  were  constrained  at  the  root  of  the  wing  in  all  six  degrees  of  freedom. 
Table  2  compares  the  models  first  five  natural  frequencies  of  the  models. 


Table  2:  Natural  Frequency  Comparison  of  Modified  and  Original  RLV  FEA  Wing 


Mode 

Original  (Hz) 

Modified  (Hz) 

%  Difference  From  Original 

1 

7.79 

8.02 

2.95 

2 

16.89 

12.57 

25.5 

3 

23.63 

26.16 

10.70 

4 

57.38 

54.88 

4.35 

5 

76.15 

72.27 

5.08 
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The  first  five  modes  were  selected  to  compare  because  in  most  cases  the  natural 
frequencies  will  most  likely  never  reach  above  them  in  the  physical  world.  Figure  6  illustrates 
the  first  five  mode  shapes  and  natural  frequencies,  where  the  figures  on  the  left  are  the  original 
wing  configuration  and  the  figures  on  the  right  are  of  the  modified  wing  configuration.  The  first 
mode  shape  of  the  wing  structures  contains  a  bending  mode.  Whereas  the  second  mode  shapes 
illustrate  a  twisting  modes.  The  third  mode  is  a  bending/twisting  mode.  The  fourth  and  fifth 
modes  are  a  combination  of  different  mode  shapes.  Although  the  natural  frequencies  do  not 
compare  exactly  some  disagreement  will  occur  when  a  lower  fidelity  model  is  introduced,  due  to 
the  reduction  of  nodes. 
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Orignal  FEA  model 


Modified  FEA  model 


MODE  1 


MODE  2 


MODE  3 


MODE  4 


MODE  5 


Figure  6:  First  Five  Natural  Frequencies  and  Mode  Shapes  of  the  Original  and  Modified  RLV  FEA  Wing  Model 
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4.2.2  Modal  Assurance  Criteria  Analysis 


Another  important  characteristic  when  comparing  two  models  is  to  identify  if  the  mode 
shapes  demonstrate  the  same  behavior.  The  modal  assurance  criterion  (MAC)  (Equation  (1))  was 
the  second  method  used  to  validate  the  new  model.  The  MAC  is  a  measure  between  two  modes, 
whether  they  are  by  two  experimental  modes,  a  physic  derived  mode  and  experimental  modes,  or 
two  physic  derived  modes  [1 1,  12].  In  most  cases  the  MAC  is  used  to  compare  a  physical 
derived  mode  and  an  experimental  mode.  Since  there  is  not  a  physical  model  of  the  RLV  to  get 
experimental  modes  the  two  finite  element  models  were  compared  assuming  the  original  model 
is  exact. 


MAC(v|/!,  v|ri) 


fil'P'I'm)2 

(4iH+l)((4m)%im) 


(1) 


The  MAC  matrix  can  be  any  size  depending  on  how  many  mode  shapes  are  being 
evaluated.  In  the  MAC  matrix  each  row  and  each  column  represent  one  mode  shape.  The  MAC 
matrix  consists  of  coefficients  between  1  and  0,  where  a  value  of  1  indicates  a  perfect  match  and 
0  indicated  no  similarity.  If  the  two  models  were  an  exact  match  there  would  be  Is  down  the 
diagonal  of  the  MAC.  A  value  of  .9  indicates  a  strong  similarity  between  the  two  shapes.  Mode 
switching  can  also  be  determined  by  the  MAC  matrix  this  can  be  seen  by  a  high  value  located  in 
an  off  diagonal. 


The  MAC  was  used  to  compare  the  first  five  mode  shapes  of  the  original  and  modified 
model.  All  six  degrees  of  freedom  of  the  mode  shapes  were  compared.  Table  3  contains  the 
results  of  the  MAC  analysis. 
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Table  3:  MAC  of  Original  and  Modified  Models 


Original  Model 

Modified 

Model 

Modes 

1 

2 

3 

4 

5 

1 

0.095065 

0.256329 

0.058633 

3.80E-05 

2 

0.057177 

0.21611 

0.069353 

0.000554 

3 

0.240514 

0.002087 

0.033638 

0.368551 

4 

0.001051 

0.140134 

0.002908 

7.78E-06 

0.737396 

5 

0.034796 

0.003098 

0.002735 

0.802372 

0.013169 

The  table  was  color-coded  to  help  distinguish  the  MAC  coefficients.  The  boxes 
highlighted  in  red  illustrate  a  good  similarity  between  the  two  shapes.  The  first  and  second  mode 
shape  of  the  models  compare  very  nicely,  which  can  be  seen  by  the  values  above  .9  on  the  table. 
This  is  critical  because  when  a  flutter  analysis  is  performed  the  first  and  second  mode  shapes  are 
usually  responsible  for  the  flutter  excitation  (which  will  be  described  in  a  later  section).  The  third 
mode  shape  is  slightly  less  accurate  but  still  correlates  with  the  original  model  demonstrated  by 
the  .87  value  in  the  diagonal.  As  described  before  in  the  previous  section  the  third  mode  shape 
displays  a  bending/twisting  motion.  This  can  be  seen  in  the  blue  boxes  that  there  is  around 
twenty  percent  correlation  between  the  first,  second,  and  third  mode  shapes.  The  fourth  and  fifth 
mode  shape  flip  and  can  be  seen  with  the  high  values  found  in  the  off  diagonal  located  in  the 
yellow  boxes. 

The  two  forms  of  validation  that  were  investigated  provides  confidence,  that  the  new 
model  developed  is  an  acceptable  surrogate  for  the  configuration  F  wing.  It  captures  the  physical 
model  closely  enough  that  the  uncertainty  quantification  analysis  can  be  implemented  with 
meaningful  data.  The  new  model  saves  computational  time  while  the  uncertainty  quantification 


20 


methods  are  being  used  and  will  allow  a  higher  fidelity  model  to  easily  replace  the  model  at  the 
appropriate  time  of  the  design  process. 

4.3  Aeroelastic  Flutter  Model  Development 

Air  vehicles  that  have  vertical  tails  on  the  tips  of  the  wing  like  the  current  RLV 
configuration  are  more  susceptible  to  flutter  due  to  the  additional  mass  on  the  wing  tip  (which 
lowers  bending  frequency)  and  the  additional  aerodynamic  lifting  surface  placed  at  the  end  of  the 
wing.  There  have  been  many  types  of  aeroelastic  analysis  conducted  on  hypersonic  and  reusable 
launch  vehicles  [13-  15].  This  provides  an  excellent  analysis  to  complete  a  risk  quantification 
study  using  aeroelastic  flutter  as  the  failure  bound  (limit  state). 

To  conduct  an  aerodynamic  flutter  analysis  many  inputs  are  needed  for  a  simulation.  The 
first  input  needed  is  a  structural  model  that  represents  the  vehicle  being  analyzed.  Next,  an 
aerodynamic  model  must  also  be  constructed  and  positioned,  to  capture  the  aerodynamic  loads 
associated  with  flight  found  in  the  trajectory  optimization.  In  this  section  the  steps  taken  to 
complete  a  series  of  flutter  analysis  on  the  RLV  wing  was  conducted. 

Aeroelasticity  is  a  primary  example  of  multi-disciplinary  analysis.  The  coupling  of 
aerodynamic  and  structural  forces  and  responses  are  analyzed  and  the  effects  that  this  has  on 
structural  response,  vehicle  perfonnance,  and  vehicle  controllability  are  able  to  be  determined. 
Perfonnance  parameters  such  as  divergent  dynamic  pressures,  control  surface  reversal  pressures, 
trim  conditions  and  flutter  velocities  can  be  calculated.  Although  aeroelastic  analyses  differ 
between  Mach  regimes  the  inputs  and  results  to  the  analyses  remain  nearly  identical  for 
implementation  within  a  design  framework. 
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Flutter  is  a  self-feeding  vibration  where  aerodynamic  forces  on  a  wing  couples  with  the 
natural  mode  of  vibration  to  produce  periodic  motion  that  can  cause  catastrophic  failure.  The 
vibration  movement  of  the  wing  increases  an  aerodynamic  load  which  in  turn  drives  the  wing  to 
deflect  further.  If  the  energy  during  the  period  of  aerodynamic  excitation  is  larger  than  the 
natural  damping  of  the  system,  the  level  of  vibration  will  increase.  The  vibration  levels  can  thus 
build  up  and  are  only  limited  when  the  aerodynamic  or  mechanical  damping  of  the  object  match 
the  energy  input,  this  often  results  in  large  amplitudes  and  can  lead  to  rapid  failure.  This  leads  to 
an  excellent  parameter  to  examine  with  the  tools  for  uncertainty  quantification  since  there  are  so 
many  uncertain  variables  that  that  could  affect  flutter. 

4.3.1  Reusable  Launch  Vehicle  Aeroelastic  Flutter  Model 

NASTRAN  FEA  software  was  selected  to  construct  the  aeroelastic  model  due  to  the 
capability  to  perfonn  linear  aeroelastic  analysis.  NASTRAN’s  aerodynamic  analysis,  like 
structural  analysis,  is  based  on  the  finite  element  approach  [16].  The  finite  aerodynamic  elements 
are  strips  or  boxes  on  which  there  are  aerodynamic  forces.  These  can  be  described  simply  by 
defining  properties  of  the  array  (panel).  The  grid  points  defining  the  structure  usually  do  not 
coincide  with  the  grid  points  defining  the  aerodynamic  elements,  provision  has  been  made  to 
generate  equations  for  interpolating  between  the  two,  known  as  splines.  This  interpolation  is  a 
key  feature,  since  it  allows  the  choice  of  structural  and  aerodynamic  considerations  to  occur 
independently.  NASTRAN  provides  three  methods  to  complete  a  flutter  analysis  including  the  K, 
P,  and  PK  methods.  All  of  the  flutter  analysis  presented  in  this  research  utilized  the  PK  method. 
Choosing  an  appropriate  aerodynamic  theory  is  critical  to  compute  a  flutter  analysis.  NASTRAN 
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has  implemented  six  aerodynamic  theories  including  doublet-lattice  subsonic  lifting  surface 
theory,  ZONA5 1  supersonic  lifting  surface  theory,  subsonic  wing-body  interference  theory, 
mach  box  method,  strip  theory,  and  piston  theory.  Doublet  lattice  and  ZONA5 1  were  both  used 
in  the  flutter  analyses  given  that  only  the  subsonic  and  supersonic  regimes  were  being  analyzed. 
The  Doublet-Lattice  method  can  be  used  for  interfering  lifting  surfaces  in  subsonic  flow. 
ZONA5 1  is  a  supersonic  lifting  surface  theory  that  accounts  for  the  interference  among  multiple 
lifting  surfaces.  Figure  7  illustrates  the  structural  and  aerodynamic  model  developed  in 
NASTRAN  for  the  aeroelastic  flutter  analysis.  The  figure  shows  the  aerodynamic  model  is  the 
same  shape  as  the  structural  model  allowing  for  accurate  interpolation  between  the  two  types  of 
models. 


Figure  7:  NASTRAN’s  RLV  Wing  Aerodynamic  Model  (Top)  and  Structural  Model  (Bottom) 


23 


4.4  FAST  Aeroelastic  Flutter  Analysis 


Figure  8  and  Figure  9  illustrate  the  results  for  a  single  aeroelastic  flutter  case,  in 
NASTRAN,  of  the  modified  RLV  wing  model,  using  Mach  1.1,  at  75  deg  F,  and  symmetric 
boundary  conditions.  The  o  represents  the  first  bending  mode,  the  □  is  the  second  mode  of  the 
wing  in  torsion,  and  the  0  is  the  third  mode  which  consists  of  a  bending  twisting  mode.  Flutter 
occurs  when  one  of  the  modes  cross  the  zero  axis  in  the  V-G  plot.  The  V-G  plot  in  Figure  8 
demonstrates  the  RLV’s  wing  begins  to  flutter  at  1500  ft/s.  The  two  modes  accountable  for  the 
flutter  are  the  first  wing  bending  and  first  torsion  modes.  Figure  9  reveals  the  coalescent  behavior 
of  the  two  modes  to  be  around  1500  ft/s  to  2000  ft/s,  also  confirming  the  flutter  is  occurring. 
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Figure  8:  RLV’s  Wing  Flutter  V-G  Plot  Representative  (Mach  1.1,  No  Fuel,  Sym.  BC,  75°  F) 


Figure  9:  RLV’s  Wing  Flutter  V-to  Plot  Representative  (Mach  1.1,  No  Fuel,  Sym.  BC,  75°  F) 
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To  further  explore  the  aeroelastic  flutter  of  the  RLV  wing  configuration  a  more  in  depth 
analysis  was  performed  to  make  an  assessment  of  flutter  speed  across  the  flight  envelope  using  linear 
aerodynamic  methods.  This  analysis  will  only  be  valid  in  the  subsonic  and  supersonic  regime  due  to 
the  panel  methods  not  being  valid  at  transonic  Mach  numbers.  This  analysis  will  use  the  same  initial 
conditions  as  stated  in  the  previous  aeroelastic  flutter  study,  only  changing  the  Mach  number. 

Figure  10  shows  the  results  for  the  full  flutter  analysis.  The  lowest  flutter  dynamic 
pressure  Mach  number  occurs  at  Mach  1.1.  In  this  analysis  the  first  mode  shape  was  responsible 
for  the  flutter  excitation  in  each  of  the  simulations. 
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5.0  RISK  QUANTIFICATION 


Uncertainty  in  a  design  problem  can  arise  from  multiple  sources  as  stated  by  Rudiger 
[17].  This  is  exacerbated  in  problems  that  involve  modeling,  whether  that  is  a  simplified 
physical  model,  or  a  computational  model.  These  uncertainties  can  commonly  associate  with 
three  distinct  sources:  Model-form  uncertainty,  parametric  uncertainty,  and  predictive 
uncertainty.  (Figure  11) 
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Figure  11:  Uncertainty  Breakdown  in  Modeling  Problems 
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The  significance  of  the  uncertainty,  and  its  source,  is  usually  related  to  the  confidence 
that  a  designer  has  in  the  models  being  utilized  for  analysis.  As  models  are  better  understood, 
such  as  with  static  finite  element  analysis,  the  contribution  to  the  overall  uncertainty  of  both 
model-form  uncertainty  and  predictive  uncertainty  are  decreased,  sometimes  to  the  point  of 
insignificance.  However,  in  less  understood  phenomenon,  or  when  working  with  over-simplified 
or  novel  models,  the  contribution  to  the  total  uncertainty  of  these  two  sources  can  potentially 
outweigh  that  of  predictive  uncertainty  by  orders  of  magnitude. 

Parametric  uncertainty  refers  to  the  uncertainty  in  the  inputs  to  the  model  and  analysis. 
Although  models  will  usually  operate  with  deterministic  parameters  in  the  scope  of  the  actual 
physical  manifestation  of  the  problem,  these  values  cannot  always  be  considered  deterministic. 
Inherent  uncertainties  can  come  from  multiple  sources,  which  contribute  and  compound  to 
produce  a  degree  of  disbelief  in  the  result  of  a  single  analysis.  While  uncertainties  in  the  design 
variables  and  input  parameters  to  a  model  are  referred  to  in  general  as  parametric  uncertainty,  the 
way  in  which  these  uncertainties  are  modeled  and  considered  determines  their  further 
classification.  Parametric  uncertainty  is  expanded  into  aleatory  and  epistemic  uncertainties  [18]. 

When  a  physical  problem  is  represented  as  a  model,  uncertainties  are  inherent  to  the 
modeling  process.  In  aeroelastic  design  these  uncertainties  can  arise  from  multiple  sources,  such 
as  the  input  parameters  into  the  model,  the  fidelity  of  aerodynamic  analyses,  and  the 
aerodynamic  and  structural  discretization  of  the  physical  domain.  Extensive  work  has  been 
completed  in  the  past  on  parametric  uncertainty  on  structural  inputs  [19,  20],  aerodynamic  inputs 
[20,  21],  and  environmental  loading  conditions  [22].  However,  the  full  uncertainty  associated 
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with  modeling  consists  of  more  than  parametric  input,  but  also  includes  other  factors,  such  as  the 
uncertainties  introduced  by  the  modeling  process  itself.  The  total  uncertainty  resulting  from 
modeling  can  at  a  high  level  be  broken  into  three  distinct  components  [23]:  model-form 
uncertainty — the  degree  of  uncertainty  between  multiple  models  of  the  same  physical  problem — 
parametric  uncertainty — the  uncertainty  in  the  parameters  of  an  analysis — and  predictive 
uncertainty — the  unknown  errors  introduced  by  the  simplifying  assumptions  of  a  model,  often 
reduced  to  be  the  difference  between  a  simulation  result  and  the  true  physical  problem  (Equation 
(2)). 


y  =  fi(x)  +  e  (2) 

In  Equation  (2),  f(x)  represents  both  the  model-form  and  parametric  uncertainty  in  the 
problem  and  e  represents  the  predictive  uncertainty.  Any  uncertainty  in  the  input  vector  x  is 
considered  parametric  uncertainty.  While  this  uncertainty  is  propagated  through  the  model,  it  is 
separable  from  model-form  uncertainty  in  a  well  understood  problem.  The  function  ft  represents 
a  model  of  the  system.  When  multiple  models  are  considered,  the  difference  between  the  values 
of  f{x)  is  considered  model-form  uncertainty.  Finally,  the  difference  between  the  model’s 
representation  of  the  system,  f{x)  ,  and  the  “true”  value  of  the  analysis,  y ,  is  called  the 
predictive  uncertainty,  i.  The  predictive  uncertainty  in  a  problem  is  a  result  of  the  assumptions 
made  in  the  modeling  process.  Parametric  uncertainty  quantification  methods  and  applications 
have  been  addressed  in  depth  in  the  existing  literature  [19  -  22,  24  -  27].  However,  the  other 
two  sources  of  uncertainty— both  model-form  and  predictive  uncertainty— are  frequently  ignored 
in  aeroelastic  problems. 
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5.1  Aleatory  Uncertainty  Quantification 


Aleatory  uncertainty  is  what  is  traditionally  thought  of  when  uncertainty  is  considered  in 
a  structural  problem.  It  refers  to  the  type  of  uncertainty  where  enough  knowledge  regarding  the 
uncertainty  of  a  parameter  is  known  such  that  a  continuous  distribution  function  of  its  values  can 
be  determined  and  assumed  to  be  valid.  Aleatory  uncertainty  is  typically  defined  in  terms  of 
probabilistic  distributions  (Figure  12). 


Figure  12:  Probability  Density  Function  in  2D 


When  described  in  this  manner,  a  probabilistic  distribution  of  output  parameters  can  often 
be  calculated,  to  which  a  probability  of  failure  can  be  prescribed  (Figure  13). 
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Figure  13:  Probability  of  Failure  Representation 


However,  all  variables  cannot  always  be  assumed  to  be  probabilistic.  While  with  well 
understood  parameters  with  an  abundance  of  data  available,  it  is  possible  to  assign  a  distribution 
to  the  parameter,  with  other  less  understood  parameters,  this  might  not  be  possible,  or  might  even 
result  in  erroneous  results.  If  a  probability  is  assigned  to  a  value  to  which  only  limited  data  is 
available,  it  is  possible  to  improperly  define  the  variable,  making  the  results  of  any  analysis 
using  that  distribution,  incorrect.  Instead  of  assigning  a  distribution  to  these  parameters  in  this 
case,  another  method  must  be  explored  for  analyzing  the  uncertainty  in  these  problems. 


Aleatory  uncertainty  must  be  accounted  for  in  the  design  process,  while  still  maintaining 
computational  efficiency.  To  quantify  aleatory  uncertainty  quantification  the  following  methods 
were  explored. 
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Monte  Carlo  Simulations 


First  order  reliability  method(  FORM) 

■  Hasofer  and  Lind  iteration  method  (HL) 

■  FORM  with  adaptive  approximations,  TANA 
Second  order  reliability  method  (SORM) 

■  Breitung’s  formulation 

■  Tvedt’s  Formulation 

■  Koyluoglu’s  Fonnulation 


The  goal  of  the  preliminary  analysis  using  aleatory  uncertainty  quantification  techniques 
is  to  find  a  sufficient  method  of  uncertainty  quantification  for  the  vehicle  design,  with  specific 
interest  initially  in  the  aeroelastic  design  and  analysis  of  an  RLV  wing.  To  complete  an  efficient 
uncertainty  quantification  analysis  it  is  important  that  a  method  is  selected  that  gives  an  accurate 
reliability  with  a  minimal  amount  of  function  evaluations,  which  is  directly  related  to  simulation 
time.  During  the  comparison  of  the  aleatory  methods,  a  closed  form  nonlinear  equation  was  used 
to  demonstrate  the  capabilities  of  the  methods  while  not  using  computational  time  on  executing 
FEA. 


Each  method  listed  above  was  analyzed  and  compared  to  a  converged  uncertainty 
quantification  Monte  Carlo  Simulation  analysis,  which  can  be  considered  in  this  case  to  represent 
an  “exact”  solution,  to  show  the  relative  accuracy  of  the  reliability  produced  by  each  method. 
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During  the  analysis  the  function  calls  were  also  recorded  to  represent  the  number  of  aeroelastic 
simulations  that  would  be  required  when  a  closed  form  equation  will  no  longer  be  available. 

The  comparative  study  examines  the  different  types  of  aleatory  uncertainty 
quantification  methods.  With  the  results,  a  decision  on  which  method  provides  the  best  accuracy 
taking  into  consideration  of  function  calls  for  this  demonstration.  For  example  if  a  method  has 
99.9%  accuracy  and  takes  a  month  to  execute  the  analysis  compared  to  a  method  that  has  97% 
accuracy  that  takes  a  week  to  execute,  a  sacrifice  will  have  to  be  made  on  accuracy  or  time.  This 
will  provide  the  designer  with  infonnation  regarding  the  trade-offs  between  accuracy  and 
computational  time. 

Equation  (3)  was  selected  as  the  closed  form  equation  to  demonstrate  the  potential  of 
each  method. 


g(x)  —  xt4  -I-  2x24  +  x1x2  —  100  (3) 

The  closed  fonn  equation  represents  a  “black  box”  where  any  type  of  simulation  can  be 
placed.  The  closed  form  equation  has  two  important  characteristics  that  are  congruent  with  an 
aeroelastic  analysis.  First,  the  closed  fonn  equation  is  highly  nonlinear;  both  xi  and  x2  are  to  the 
fourth  power  on  behalf  of  the  nonlinearity  of  an  aeroelastic  analysis.  The  second  attribute  of  the 
equation  is  xi  and  x2  are  coupled,  representing  the  coupling  of  variables  between  aerodynamic 
and  structural  portions  of  the  aeroelastic  analysis.  Variables  xi  and  x2  both  have  a  mean  of  10 
and  a  normal  standard  deviation  of  5  throughout  the  comparison  study.  g(x)  represents  the  limit 
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state,  which  signifies  the  failure  boundary,  if  the  simulation  results  in  a  negative,  it  is  considered 
a  failure. 

5.1.1  Monte  Carlo  Aleatory  Uncertainty  Quantification 

Monte  Carlo  is  one  of  the  most  basic  sample  techniques  developed  by  Neuman  and  Ulam 
in  1949  known  as  a  simple  random  sampling  method.  This  technique  can  be  implemented  to 
solve  uncertainty  quantification  analyses.  To  use  the  Monte  Carlo  Simulation  technique,  a 
distribution  type  is  needed  for  the  random  variables.  With  the  distribution  type  a  random  sample 
set  is  compiled.  Using  the  values  in  the  sample  set  as  input  values  simulations  are  then  executed. 
In  most  cases  Monte  Carlo  simulation  provides  accurate  results  with  very  high  computational 
cost.  Using  the  random  variables  and  closed  form  equation  (Equation  (3)),  a  Monte  Carlo 
Simulation  was  performed  using  Matlab  mathematical  program  to  process  the  infonnation. 

The  number  of  simulations  was  increased  during  each  trial  in  the  analysis  to  improve 
accuracy  and  demonstrate  convergence.  During  each  execution  of  the  Monte  Carlo  simulation, 
the  failures  were  recorded.  The  probability  of  failure  was  calculated  by  taking  the  number  of 
failures  and  dividing  it  by  the  number  of  simulations  as  seen  in  Table  4. 
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Table  4:  Monte  Carlo  Simulations  Results 


No.  of  Simulations 

1,000 

10,000 

50,000 

100,000 

200,000 

500,000 

1,000,000 

2,000,000 

3,861,000 

Number  of  failures 

5 

48 

234 

468 

956 

2379 

4769 

9476 

18335 

P 

f 

(probability  of  failure) 

0.005 

0.0048 

0.00468 

0.00468 

0.00478 

0.004758 

0.004769 

0.004738 

0.00474 

Noting: 

1 .  The  number  of  simulations  correlates  to  the  accuracy  of  the  Monte  Carlo  Simulation.  As 
seen  in  Table  4  the  probability  of  failure  converges  as  more  simulations  are  conducted. 

2.  A  convenience  of  the  Monte  Carlo  sampling  technique  is  the  process  can  be  stopped  and 
started  at  any  time  because  each  simulation  is  independent  from  the  next. 

3.  Knowing  each  simulation  is  independent  from  the  last,  saved  time  in  this  comparative 
study. 


The  final  simulation  of  the  Monte  Carlo  technique  employs  3,861,000  function 
evaluations.  The  analysis  resulted  in  the  probability  of  failure  of  0.00474,  and  will  be  used  as  the 
most  accurate  uncertainty. 


5.1.2  First  Order  Reliability  Methods  (FORM) 

As  seen  in  the  results  from  the  Monte  Carlo  Simulation,  many  simulations  are  needed  to 
achieve  a  respectable  solution;  this  provides  motivation  to  find  a  more  efficient  method.  First 
order  reliability  method  (FORM)  is  an  alternative  technique  to  solve  aleatory  uncertainty 
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problems.  FORM  uses  an  iteration  method  while  approximating  the  limit  state  function  to  solve 
for  the  probability  of  failure.  Yan-Gang  et  al  [28]  demonstrated  a  general  procedure  for  the 
first/second-order  reliability  methods  will  be  used  in  this  investigation. 

Figure  14  graphically  represents  the  FORM  process.  The  first  step  in  the  process  is  the 
same  as  Monte  Carlo,  where  the  mean  values,  standard  deviations,  and  distributions  of  the 
uncertain  variables  are  required.  In  the  figure  pi  and  po  both  have  normal  distributions.  g(x)  =  0 
represents  the  limit  state  which  separates  the  safe  and  failure  region.  FORM  requires  the  means 
of  the  uncertain  variables  and  the  limit  state  function  be  normalized  in  this  case  using  Equation 
(4)  (X-space  is  transfonned  to  U-space). 
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Figure  14:  Graphical  FORM  Representation 


Ui  = 


Xi  M  X ; 


(4) 


37 


This  transformation  allows  for  the  means  of  the  uncertain  variables  to  relocate  to  the 


origin  (0,  0  position).  Once  transformed  to  U-space  an  optimization  is  preformed  to  minimize  the 
length  of  (3  (reliability  index).  Where  (3  contacts  the  limit  state  function  the  most  probable  point 
of  failure  is  detennined  (MPP).  A  first  order  approximation  can  be  made  at  the  MPP.  This 
approximation  is  illustrated  by  a  line  drawn  tangent  to  the  limit  state  function  at  the  MPP.  The 
reliability  is  then  found  by  locating  the  correlating  percent  in  a  cumulative  standard  normal 
distribution  table  to  determine  the  probability  of  failure.  Two  methods  to  solve  FORM  problems 
are  Hasofer  and  Lind  iteration  method  (HL)  and  FORM  with  adaptive  approximations. 


The  Hasofer  and  Lind  iteration  method  (HL)  is  a  recursive  algorithm  used  to  solve 
reliability  problems  [29].  The  HL  method  approximated  the  limit-state  function  using  the  first- 
order  Taylor  series. 

FORM  with  adaptive  approximation  also  approximates  the  limit  state  function  using  two- 
point  adaptive  nonlinearity  approximation  (TANA)  Equation  (5).  Although  extra  function  calls 
are  needed  initially  to  use  TANA  it  closely  approximates  the  limit  state  function  resulting  in  less 
function  calls  to  convergence  [30], 

GTm  =  G(u2 )  +  -  u,/)  g-)  (5) 

l'  U 2 
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5.1.3  Second  Order  Reliability  Methods  (SORM) 


FORM  works  well  when  the  limit  state  function  is  nearly  linear.  However,  when 
considering  aeroelastic  limit  state  the  function  is  not  always  linear.  A  second  order  reliability 
method  (SORM)  will  tend  to  provide  a  more  accurate  representation  of  the  design  space.  SORM 
requires  the  p  and  MPP  from  the  FORM  analysis.  SORM  calls  for  second  order  derivatives  of 
the  limit  state  function.  On  top  of  Figure  15  a  first  order  approximation  is  shown.  As  observed  in 
the  graph  the  first  order  approximation  underestimates  the  probability  of  failure  by  a  significant 
amount,  giving  inaccurate  results.  On  the  bottom  of  Figure  15  a  second  order  approximation  is 
used  to  illustrate  the  improvement  of  the  probability  of  failure  approximation. 


First-order  Approximation  ofthe  response 
Surface  in  Y-Space 


Limit  state  function 


Second-order  Approximation  ofthe  response 
Surface  in  Y-Space 

Figure  15:  Graphical  Representation  of  FORM  vs.  SORM 
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The  first  step  of  the  process  is  to  find  P  and  the  MPP  in  U  space  using  a  FORM  method. 
The  next  step  is  to  compute  the  second-order  derivatives  and  create  the  B  matrix  which  becomes 
very  computationally  expensive  when  more  than  two  variables  are  being  examined.  The 
computational  time  can  be  reduced  by  using  a  function  approximation  (Such  as  TANA  described 
in  the  previous  section).  Next  the  H  matrix  is  created  by  rotating  U  space  to  Y  space.  This  is 
completed  by  orthogonalizing  the  H  matrix  using  the  Gram  Schmidt  algorithm.  Following  with 
computing  the  main  curvatures  k,  by  solving  the  eigenvalues  of  HBHT.  k,  and  P  is  then  used  in 
each  of  the  following  methods  to  compute  the  probability  of  failure.  The  following  SORM 
methods  have  a  brief  explanation  and  show  the  equations  for  the  probability  of  failure. 

Breitung’s  formulation  was  one  of  the  first  methods  introduced  in  1984  [31],  Equation 
(6).  The  formulation  is  derived  as  an  asymptotic  formula  of  the  failure  probability. 


Vf  =  *(-/?/)  npil+Pfkj)1  (6) 

Tvedt  introduced  a  three-tenn  approximation  in  1984  in  which  A2  and  A3  can  be 
interpreted  as  the  correction  for  Breitung’s  formula  [32]  Equation  (7). 
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(7) 


Pf  -  'K-/?/)  riy=i  (l  +  Pfkj)  2  +A2  +A3 


a2  =  -  </>(-%)]  {n;=i(i + mP  -  +  wf  +  m  pj 

A3  =  (Pf  +  +  WP  -  He  nj'rJCl  +  (Pf  +  1)*;P} 


Koyluoglu  and  Nielsen  developed  Koyluoglu’s  formulation  in  1994  making  it  a  more 
recent  method  compared  to  the  others  [33]  Equation  (8).  This  fonnulation  is  a  one-tenn 
approximation. 


71  —  1 


for  k:  >  0 ,j  —  1,2,  ...,71-1 


for  ki  <  0 ,j  —  1,2, ... ,  n  —  1 


(8) 


Using  FORM  and  SORM  methods,  an  uncertainty  quantification  analysis  was 
accomplished  on  the  closed  fonn  equation  represented  by  Equation  (3)  using  the  same  conditions 
that  were  used  in  the  Monte  Carlo  simulation.  The  results  of  this  comparative  study  can  be  found 
in  Table  5.  The  gradients  of  the  closed  form  equation  are  found  using  finite  difference  rather  than 
mathematically  solving  for  them  to  represent  a  black  box  simulation. 
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Table  5:  UQ  Results  for  a  Closed-Form  Problem 


Total  Number  of  Function 

Calls 

Pf 

%  Difference  from 

MC 

HL  iteration  method 

87 

0.01589 

235.4 

FORM  with  adaptive 
approximations 

15 

0.01590 

235.4 

Breitung's  Formulation 

21 

0.00498 

5.0 

Tvedt's  Fonnulation 

21 

0.00436 

7.95 

Koyluoglu's  Formulation 

21 

0.00464 

2.08 

Monte  Carlo 

3,861,000 

0.00474 

0 

The  results  from  the  aleatory  uncertainty  quantification  methods  were  compared.  The  HL 
iteration  method  and  FORM  with  adaptive  approximation  perfonned  poorly  as  expected  because 
the  limit  state  function  was  nonlinear.  FORM  with  adaptive  approximations  using  TANA  as  the 
approximation  shows  a  much  faster  convergence  than  the  HL  iteration  method.  The 
improvement  of  FORM  with  adaptive  approximations  finds  an  optimum  p  in  15  functions  which 
is  directly  related  to  the  amount  of  computation  time  that  it  would  take  when  using  this  method 
for  vehicle  design.  The  FORM  with  adaptive  approximations  method  will  be  used  during  the 
uncertainty  quantification  in  vehicle  design  based  on  the  savings  of  computational  time. 

Breitung's  Formulation,  Tvedt' s  Fonnulation,  and  Koyluoglu's  Fonnulation  were  the 
three  SORM  methods  compared.  As  explained  earlier  a  FORM  analysis  is  needed  to  construct 
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SORM.  This  makes  SORM  have  a  higher  computational  cost  then  FORM.  Since  the  same 
method  is  used  to  solve  for  the  Eigenvalue  k  that  is  used  in  the  formulas  the  three  methods  have 
the  same  amount  of  function  calls.  For  this  particular  problem  Koyluoglu's  Formulation 
performed  closest  to  the  Monte  Carlo  results. 

5.2  Epistemic  Uncertainty  Quantification 

In  problems  where  little  variable  infonnation  exists,  it  is  necessary  to  use  an  uncertainty 
quantification  method  that  does  not  introduce  additional  assumptions,  which  would  result  in 
adding  more  uncertainty  into  the  problem.  These  methods  are  referred  to  as  epistemic 
uncertainty  quantification  methods.  Epistemic  uncertainty  is  uncertainty  to  which  no 
assumptions  are  made  regarding  the  parameters.  Instead  of  defining  a  distribution  to  the 
parameters,  intervals  and  bounds  of  parameters  are  assigned  based  upon  limited  available  data, 
expert  opinions,  or  prior  knowledge  of  the  problem.  Using  these  interval  definitions  of  variables, 
the  uncertainty  in  the  system  must  then  be  propagated  using  advanced  methods.  A  known 
method  of  propagating  this  uncertainty  is  Dempster-Shafer  Theory  or  Evidence  Theory  [34,  35], 

5.2.1  Evidence  Theory 

A  technique  is  needed  to  quantify  uncertainty  when  there  is  little  infonnation  known 
about  the  uncertain  variables.  Two  epistemic  variables  that  play  a  key  role  in  a  RLV  are 
atmospheric  conditions  and  composite  properties.  Evidence  theory  can  quantify  epistemic 
uncertainty  without  making  any  assumptions  and  was  introduced  to  the  reliability  of  structures 
by  Bae  et  al  [26,  36].  Evidence  Theory  measures  uncertainty  with  two  measures,  belief  (BEL) 
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and  plausibility  (PL).  The  reliability  of  the  problem  can  be  found  between  these  two  bounds.  The 
upper  bound  is  the  plausibility  and  the  lower  is  belief,  this  defines  a  probability  interval.  The 
bound  size  changes  by  the  amount  of  infonnation  known  about  the  uncertain  variables.  When 
little  infonnation  is  known  the  bound  is  large,  contrary  when  full  information  is  known  the 
bound  is  small.  To  establish  the  values  of  BEL  and  PL  a  basic  belief  assignment  (BBA)  (Figure 
17)  needs  to  be  constructed.  To  construct  a  BBA,  first  a  frame  of  discernment  needs  to  be 
established.  The  frame  of  discernment  represents  (Figure  16)  all  the  possible  distinct 
propositions  given  in  Equation  (9). 


Xj  x2  x3 

Figure  16:  Frame  of  Discernment 


2X  =  (0,  {xt},  {x2},  {x3},  (xlf  x2j,  {x2,  x3},  {xt,  x3],  X}  (9) 

Each  combination  in  Equation  (9)  can  be  assigned  a  value  from  zero  to  one.  This  value  will  be 
referred  to  as  m  and  is  the  weight  of  belief  for  that  section  of  the  frame  of  discernment. 


m 


2 


Figure  17:  Basic  Belief  Assignment 
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The  m  contains  the  weighted  information  about  the  BBA  and  must  satisfy  Equations  (10  -  12). 


m(A)  >  Of  or  any  A  G  2X  (10) 

m(0)  =0  (11) 

ZAe2xm(A)  =  1  (12) 

Using  the  BBA  structure  defined,  Equations  (13,  14)  define  BEL(A)  and  PL(A)  respectively. 

BEL(A)=ZCcAm(Q  (13) 

PL(A)  =  ZcnA*om(C)  (14) 

To  demonstrate  evidence  theory,  a  two  variable  case  is  provided  based  on  an  RLV  wing 
flutter  reliability.  In  this  analysis  the  limit  state  is  2000  pounds  per  square  foot  (psf)  flutter 
dynamic  pressure.  This  states  that  if  the  flutter  dynamic  pressure  falls  below  2000  psf  the 
analysis  is  considered  a  failure.  Figure  18  denotes  the  two  BBA’s  for  the  two  variable 
demonstration  case.  The  uncertain  variables,  temperature  and  gas  constant,  (which  relate  to  the 
uncertainties  with  the  air  density)  are  considered  epistemic  variables  because  there  is  not  enough 
accurate  infonnation  to  create  a  valid  pdf.  In  this  case  the  intervals  were  chosen  by  experts  and 
the  m  was  assigned  by  the  expert’s  opinion. 
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Figure  18:  Gas  constant  and  Temperature  BBA  for  Evidence  theory  Demonstration  Case 


The  next  step  is  to  combine  the  variables  and  complete  the  analysis.  Figure  19  shows  how 
the  variables  interact  with  one  another.  This  is  done  by  putting  BBA’s  on  the  x  and  y  axis.  Lines 
are  then  drawn  at  each  interval  to  show  where  the  variables  intersect.  The  intersections  are 
represented  by  black  dots.  At  each  black  dot  a  flutter  analysis  was  executed  with  the 
corresponding  inputs  found  on  the  x  and  y  axis.  This  two  variable  problem  shows,  with  the 
BBA’s  provided,  20  simulations  will  need  to  be  executed.  To  calculate  the  BEL  and  PL 
hypercubes  must  be  calculated.  A  hypercube  in  a  two  variable  problem  is  represented  by  an  area. 
For,  example  the  yellow  box  in  Figure  19  is  a  hypercube  in  which  two  of  the  intervals  interact 
with  one  another.  The  area  is  then  calculated  by  multiplying  .6  by  .1  which  results  in  m=.06.  The 
hypercubes  must  be  calculated  for  each  interaction. 
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Figure  19:  Two  Variable  Hypercube  Calculation 


Once  the  hypercubes  are  calculated  the  next  step  can  be  seen  in  Figure  20.  The  Flutter 
analysis  was  executed  at  each  intersection.  If  the  flutter  analysis  resulted  in  a  flutter  dynamic 
pressure  less  than  2000  psf  the  point  was  represented  with  a  red  dot.  Likewise,  if  the  flutter 
resulted  in  a  dynamic  pressure  more  than  2000  psf  the  intersection  was  represented  with  a  blue 
dot. 
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Figure  20:  Two  Variable  Demonstration  Case  Function  Evaluation 


The  BEL  and  PL  can  be  calculated  by  examining  Figure  20.  If  all  of  the  intersections  of  a 
hypercube  are  blue  then  that  hypercube  is  part  of  the  BEL  equation.  If  a  hypercube  has  a 
combination  of  red  and  blue  dots  it  falls  in  the  PL  section.  PL  also  includes  all  of  the  hypercubes 
that  were  in  the  BEL  equation.  If  a  hypercube  contains  all  red  dots  it  is  considered  a  failure 
region. 
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BEL(A)  =  V  m(C) 

Cc.4 

BEL=.015+.03+.51 


BEL=.555 


PL(A)  =  I  m(C) 
crvi^o 

PL=.015+.03+.51+.225+.03+.06 

PL=.87 


Figure  21:  Two  Variable  Belief  and  Plausibility  Calculations 


Figure  21  shows  the  calculation  of  BEL  and  PL.  The  BEL  and  PL  are  the  bounds  of  the 
probability  of  success.  This  illustrates  without  making  any  assumptions  that  the  reliability  falls 
between  .555  and  .87. 

5.3  Plausibility  Decision 

To  accomplish  the  task  of  merging  evidence  theory  into  mathematical  optimization 
procedures,  sensitivity  infonnation  of  BEL  and  or  PL  metrics  must  be  available  with  respect  to 
the  design  variables.  Design  optimization  using  a  search  direction  is  one  of  the  most  basic 
methods  of  structural  optimization.  The  most  straightforward  approach  to  obtain  gradients,  with 
respect  to  the  design  variables,  is  the  finite  difference  method.  As  seen  in  Figure  22  a  belief  and 
plausibility  demonstration  is  plotted  with  respect  to  a  design  variable.  This  shows  even  if  finite 
difference  was  used  there  is  little  change  in  the  belief  and  plausibility  when  the  design  variable  is 
altered.  The  resulting  gradients  would  always  be  a  one  or  zero  making  gradient-based  design 
optimization  nearly  impossible. 
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Figure  22:  PL  and  BEL  Discontinues  Behavior  Based  on  Design  Variable  Modification 


To  use  gradient-based  design  optimization,  a  new  measure  must  be  detennined.  Bae 
introduced  a  new  measure  known  as  plausibility  decision  (PLdec  )  [37].  Alyanak  introduced  an 
additional  approximation  method  to  find  PL  dec  [38].  This  measure  assumes  a  uniform 
probability  distribution  for  each  distinct  proposition  interval  in  the  frame  of  discernment.  Then 
the  uncertainty  is  directly  integrated,  using  approximation  functions  to  increase  efficiency.  Since 
PL  dec  is  detennined  by  integration,  it  creates  a  continuous  measure  that  can  be  seen  in  Figure 
23.  Given  that  PL  dec  is  a  continuous  measure  it  will  be  possible  to  find  gradient  infonnation 
that  can  be  used  in  a  finite  difference  method. 
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Figure  23:  Plausibility  Decision  Behavior  Based  on  Design  Variable  Modification 


Three  methods  presenting  the  plausibility  decision  approximation  can  be  found  in  Figure 
24.  On  the  left  is  the  approximation  developed  by  Bae  in  which  the  limit  state  function  is 
approximated  by  using  a  nonlinear  approximation  (TANA).  In  the  middle  Alyanak  developed  a 
linear  approximation  that  reduces  computational  time.  Benanzer  developed  a  numerical 
approximation  where  no  extra  simulations  are  required  [39].  Benanzer  plausibility  decision 
method  is  an  internal  method  used  at  Wright  State  University  (Equation  (15)). 


Bae  [37]  Alyanak  [38]  Benanzer  [39] 

Figure  24:  Plausibility  Decision  Methods 
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X,V  !  max -(c/,,0) 

YtiJa? 


(15) 


Each  of  the  methods  has  previously  been  validated  and  can  be  implemented  into  a  finite 
difference  scheme  to  find  the  gradients  of  uncertainty  quantification  analysis.  One  major 
difference  in  the  approximations  is  how  they  behave  in  an  n-dimensional  problem.  The 
complexity  of  the  approximation  of  PL  dec  exponentially  expands  when  more  variables  are  used 
in  Bae’s  and  Alyanak’s  approximations.  One  advantage  to  Benanzer’s  plausibility  decision 
approximation  is  it  can  easy  be  implemented  into  an  n-dimensional  problem,  because  the 
approximation  is  completed  numerically.  A  validation  of  Benanzer’s  plausibility  decision 
approximation  using  two  uncertain  variables  in  one  hypercube  that  contains  plausibility,  can  be 
seen  in  Figure  25  using  Equation  (16)  as  the  limit  state.  The  results  of  the  actual  and 
approximation  plausibility  decision  can  be  found  in  Table  6.  From  this  point  onward  Benanzer’s 
plausibility  decision  approximation  will  be  used  in  all  of  the  uncertainty  quantification  analyses. 


(16) 
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Figure  25:  Two  Variable  Hypercube  Demonstrating  Plausibility  Decision 


52 


Table  6:  Two  Variable  Plausibility  Decision  Results 


Simulation 

xl_l 

xl_2 

x2_l 

x2_2 

Actual 

PL_Dec 

%  Difference 

1 

9.3 

9.7 

9.3 

9.8 

0.2947 

0.3155 

7.05 

2 

9.2 

9.5 

9.1 

9.5 

0.9769 

0.9425 

3.52 

3 

8.9 

9.6 

8.7 

9.8 

0.8797 

0.8512 

3.24 

4 

8.7 

10.3 

8.8 

10.8 

0.3801 

0.4374 

15.07 

5 

8.3 

9.3 

8.8 

11.8 

0.7005 

0.7371 

5.22 

6 

8.8 

9.1 

10.1 

11.7 

0.26 

0.2891 

11.19 

The  comparison  of  the  approximation  and  exact  plausibility  decision  is  located  in  Table 
6.  The  assessment  shows  that  the  actual  and  PLdec  range  from  three  to  fifteen  percent 
difference.  Although  this  may  seem  high  plausibility  decision  is  calculated  in  each  hypercube 
that  contains  plausibility  meaning  this  approximation  can  over  estimate  in  some  and 
underestimate  in  others  making  them  a  close  approximation  for  the  entire  evidence  theory 
analysis.  PL  dec  is  only  an  approximation  and  is  used  only  to  get  a  feel  if  the  reliability  of  a 
system  is  close  to  the  plausibility  or  belief  bound. 


To  further  expand  the  validation  of  Benanzer’s  plausibility  decision  approximation  three 
variables  in  one  hypercube  that  contains  plausibility  will  be  validates  and  can  be  seen  in  Figure 
26  using  Equation  (17)  as  the  limit  state. 


g(xllx2)  =  3 


18  _  6V3  _  2 
Xi  x2  *3 


(17) 
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Table  7  contains  the  results  from  the  three  variable  demonstration  case.  Again  this 
method  produces  a  reasonable  approximation.  In  simulation  six  it  shows  a  twenty  five  percent 
difference  of  the  approximation.  Taking  a  closer  look  .1136  and  .  1428  are  being  compared.  The 
approximation  is  only  .0292  off  of  the  actual,  but  since  the  values  are  so  low  the  percent 
difference  results  in  a  large  difference  that  is  misleading. 


Table  7:  Three  Variable  Plausibility  Decision  Results 


Simulation 

xl l 

xl 2 

x2 l 

x2 2 

x3 l 

x3 2 

Actual 

PL Dec 

%  Difference 

1 

9.1 

9.7 

9.2 

9.9 

25.1 

25.9 

0.9825 

0.9566 

2.64 

2 

9.3 

9.9 

9.2 

9.9 

22.3 

25.4 

0.8606 

0.8397 

2.43 

3 

9.2 

10.1 

9.2 

9.7 

21.3 

31.4 

0.7461 

0.7745 

3.81 

4 

9.7 

12.1 

7.4 

8 

8.4 

9.7 

0.9999 

0.9964 

0.35 

5 

10.3 

13.1 

7.9 

9.3 

7.4 

10.2 

0.4499 

0.528 

17.36 

6 

12.7 

13.1 

7.6 

8.3 

7.6 

8.5 

0.1136 

0.1428 

25.70 
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5.3.1  Sensitivity  Analysis  in  Evidence  Theory 


To  demonstrate  a  sensitivity  analysis  using  evidence  theory  and  PL  dec,  a  four  bar-truss 
structure  will  be  investigated  (Figure  27).  The  four  bar  truss  structure  was  developed  by  Haftka 
et  al  [40]  to  demonstrate  a  deterministic  constrained  optimization  problem.  The  original  problem 
has  been  modified  to  demonstrate  plausibility  decision  and  later  reliability  design  optimization. 
There  are  three  bars  in  the  structure  with  the  same  length  and  cross  sectional  area  Ai  and  a  final 
bar  with  cross  sectional  area  A2.  The  modulus  of  elasticity  is  denoted  by  E,  the  load  is  P,  and  L 
represents  the  length  of  the  beam.  The  constraint  considered  in  this  analysis  is  the  displacement 

constraint  situated  at  Equation  (18)  where  g  (x)  =  yy  (Jy  +  ■  When  the  problem  is 

converted  to  a  non-dimensional  problem  the  limit  state  can  be  seen  in  Equation  (19)  where 
x1  —  10“3  andx2  —  10-3 

I  P  z  p 
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8  <  10  ~3L 


(18) 


9  (xi.x2) 


18  _  6V3 
x1  x2 


(19) 


Using  the  BBAs  seen  in  Figure  28  for  the  variables  Xi  and  X2  a  reliability  analysis  was 
conducted.  Notice  the  BBAs  are  now  in  percentages.  The  percentages  of  the  BBA  are  needed  to 
perfonn  the  finite  difference  gradient  evaluation. 


Xi 


95%  96%  97%  98%  99%  100%  101%  102%  103% 

< - X - X - > 


.10 


.35 


.40 


.15 


x2 


95%  96%  97%  98%  99%  100%  101%  102%  103% 

< - X - X — X - > 

.80  .03  .07  .10 


Figure  28:  Four-Bar  Truss  Structure  BBAs 


To  demonstrate  the  gradients  of  belief,  plausibility,  and  plausibility  decision  with  respect 
to  the  design  variables  of  the  displacement  limit  state,  two  parametric  investigations  were 
completed.  This  is  completed  by  holding  one  of  the  design  variables  constant  at  the  design 
condition  and  sweep  the  other  design  variable  to  determine  the  exact  behavior  of  the  measures 
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Success  Estimation 


with  respect  to  the  design  variable  changes.  The  results  can  be  found  in  Figure  29  and  Figure  30 
where  the  red  line  is  the  BL,  the  blue  line  is  PL  and  the  green  dashed  line  is  PL  dec. 
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The  figures  from  the  parametric  study  show  that  belief  and  plausibility  are  discontinuous. 
The  plausibility  decision  is  a  smooth  and  continues  curve  that  falls  between  the  belief  and 
plausibility  and  gradient  information  can  be  used. 

Now  it  is  possible  to  complete  gradient-based  sensitivities  using  finite  difference  and 
plausibility  decision.  Equation  (20)  represents  the  gradient  equation. 


dP  Idee 
dxt 


(20) 


Forward  finite  difference  will  be  demonstrated  on  the  four-bar  truss  structure.  The  three 
points  ran  in  evidence  theory  using  the  four-bar  truss  structure  can  be  found  in  Table  8. 
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Table  8:  Finite  Difference  Numerical  Information 


Point 

BEL 

PL 

PLdec 

(9.5,  9.5) 

.0550 

.2775 

.1123 

(9.6,  9.5) 

.0550 

.6095 

.1975 

(9.5,  9.6) 

.0550 

.2775 

.1582 

Equations  (2 1 ,  22)  show  the  gradients  using  finite  difference  with  respect  to  the  design 
variables  Xi  and  X2.  From  the  results  it  can  be  seen  that  variable  Xi  is  almost  as  twice  as 
sensitive  as  variable  X2  at  this  design  point.  This  means  Xi  has  a  much  larger  role  in  the  design 
of  the  truss  structure  at  this  particular  point. 


dP  Ldec 
dx\ 


.1975— .1123 
9.6— 9.5 


=  .8520 


(21) 


dP  Ldec 
dxi 


.1582 -.1123 
9.6— 9.5 


=  .4590 


(22) 


5.4  Model-Form  Uncertainty  Quantification 

When  an  analysis  must  be  completed  for  an  engineering  problem,  a  representative  model 
is  often  constructed  to  allow  for  analysis  of  the  system.  To  construct  this  model,  assumptions 
regarding  the  system  must  be  made  to  simplify  the  problem  to  a  level  at  which  a  model  can 
feasibly  and  efficiently  be  constructed.  These  assumptions  often  vary  between  models  and 
modeling  packages,  resulting  in  multiple  solutions  to  identical  problems.  The  difference  between 
the  multiple  models  of  the  same  problem  is  representative  of  model-fonn  uncertainty — the 


59 


uncertainty  induced  by  the  disagreement  among  multiple  models  of  the  same  phenomenon. 
Because  there  are  multiple  models  that  give  different  answers  to  the  same  problem,  a  method 
must  be  utilized  to  combine  these  individual  results  into  a  unified  solution  while  quantifying  the 
uncertainty  associated  with  this  solution  induced  by  the  disagreement  between  the  models. 
Multiple  methods  have  been  developed  and  implemented  in  the  literature  to  quantify  this  model- 
form  uncertainty,  such  as  Bayesian  Model  Averaging  [41]  and  the  adjustment  factors  approach 
[42]. 


The  adjustment  factors  approach  was  first  demonstrated  as  a  method  to  utilize  expert 
opinions  in  Bayes’  Theorem  by  Mosleh  and  Apostolakis  [42]  in  1986.  This  method  uses  an 
adjustment  factor  to  modify  the  result  of  the  best  model.  The  adjustment  factors  approach  has 
been  demonstrated  on  multiple  engineering  problems  by  Zio  and  Apostolakis  [43]  and  Reinert 
and  Apostolakis  [44]. 

The  adjustment  factor  can  be  represented  by  multiple  types  of  distribution — such  as  a 
nonnal  or  log-normal  distribution — resulting  in  different  adjustment  factors  being  used.  In  the 
additive  adjustment  factor  approach,  the  adjustment  factor  E*a  is  assumed  to  be  a  normal  random 
variable.  The  representation  of  the  adjusted  model  prediction  is  shown  in  Equation  (23),  where 
y*  represents  the  best  model  based  on  the  expert  opinion. 

y  =  y+K  (23) 
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Knowing  the  results  of  multiple  models,  as  well  as  their  probabilities  based  upon  the 
provided  expert  opinion,  the  means  and  variances  of  the  adjusted  model  can  be  calculated  by 
Equations.  (24  -  27): 


=  (24) 

i= 1 

Var{E'a )  =  X/'fVMO-, -E(y))2  (25) 

/= 1 

E(y)  =  y+E(E'a)  (26) 


Var(y)  =  Var(E*a )  (27) 


In  the  above  equations,  £  (y)  represents  the  mean  value  of  y,  N  is  the  total  number  of 
models  considered,  P{M{)  represents  the  probability  of  model  Mt  based  upon  expert  opinion,  and 
Yi  represents  the  prediction  of  model  Mt .  From  the  above  equations,  the  mean  and  the  standard 
distribution  of  the  adjusted  model  y  can  be  calculated  and  a  normal  distribution  of  the  output  can 
be  constructed. 

In  the  multiplicative  adjustment  factors  approach,  the  adjustment  factor  E*n  is  assumed  to 
be  a  lognormal  random  variable,  and  thus  the  adjusted  model  prediction  is  shown  in  Equation 
(28). 
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7  =  7 


(28) 


Assuming  again  that  the  results  and  the  probabilities  of  the  multiple  models  are  known, 
the  means  and  variances  of  the  natural  log  of  the  adjusted  model  and  adjustment  factor  can  be 
calculated  by  Equations  (29  -  32): 


£(ln  <, )  =  Z  P(M<  XH*  I  - 1”  /  > 


i=l 


(29) 


Var(\nE*a  )  =  |>(M;)(ln|y,.|  -  £(ln|y|))2 


(30) 


E{\n\y\)  =  \ny  +E(\n  El) 


(31) 


Var(\n\y\)  =  Var(\n\El\)  (32) 

The  adjustment  factors  approach  produces  a  statistical  distribution  of  the  adjusted  model, 
accounting  for  the  variation  among  the  individual  models.  This  distribution  is  dependent  upon 
the  expert  opinion  that  goes  into  the  model  probabilities.  The  expert  opinions  are  not  infallible, 
and  thus,  an  additional  degree  of  uncertainty  is  introduced  to  the  final  distribution  due  to  the 
uncertainty  surrounding  the  model  probabilities. 
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The  uncertainty  in  the  model  probabilities  can  lead  to  multiple  problems  with  the 
adjustment  factors  approach.  Due  to  the  weighting  of  the  adjustment  factors  by  the  model 
probabilities,  it  is  possible  for  the  model  probabilities  to  have  significant  effects  upon  the 
adjusted  model.  If  changes  in  model  probabilities  lead  to  large  changes  in  the  adjusted  model, 
the  adjusted  model  becomes  very  dependent  upon  the  model  probabilities  in  addition  to  the 
variance  among  the  models.  Thus,  it  is  critical  to  be  able  to  identify  the  problems  in  which  the 
adjusted  model  is  sensitive  to  the  model  probabilities.  To  do  this,  the  Modified  Adjustment 
Factors  Approach  was  developed  in  this  work. 

The  probabilities,  assigned  to  each  model  are  initially  based  upon  expert  opinion, 

or  an  incomplete  set  of  preliminary  data.  The  uncertainty  involved  in  assigning  these 
probabilities  comprises  the  model  uncertainty  in  the  problem,  and  introduces  an  additional  layer 
of  uncertainty  in  to  the  output  distribution,  y.  To  propagate  this  uncertainty  through  to  the 
distribution  v,  the  values  P(Mt )  are  treated  as  uncertain  variables  with  a  defined  nonnal 
distribution  (Equation  (33)). 


P(Mj)  =  N(P(Mt )  exp ,  <jpm  )  (33) 


While  incorporating  the  distributions,  the  probabilities  of  each  of  the  models  must  be 
constructed  to  still  satisfy  Equation  (34). 


Ew= 1 

1=1 
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Such  that 


0  <  P(Mt )  <  1  (34) 

Due  to  the  constraints  of  Equation  (34),  one  model  probability  must  be  calculated  based 
upon  the  results  of  the  others.  Due  to  the  linear  nature  of  Equation  (34),  though,  if  the  standard 
deviations  of  all  P(Mt)  are  equal,  the  resulting  distribution  for  the  final  model  will  still  be  a 
nonnal  distribution  around  its  mean  value  with  the  same  standard  distribution  as  the  other 
models. 

As  P(Mj)  is  no  longer  a  detenninistic  measure,  its  uncertainty  must  be  propagated 
through  both  adjustment  factors  approaches  (Equations  (24  -  27),  (29  -  32)).  By  utilizing  a  Monte 
Carlo  Sampling  method,  the  statistical  data  regarding  the  final  variable  distribution  can  be 
calculated  at  individual  model  probabilities.  These  data  points  can  then  be  combined  to  form  a 
final  distribution  of  the  metric  of  interest  which  incorporates  the  uncertainty  involved  in  the 
model  probability  selection.  By  comparing  this  new  probabilistic  adjustment  factors  approach 
distribution  to  the  deterministic  distribution  calculated  by  a  traditional  approach,  a  rough 
sensitivity  of  the  global  problem  to  model  uncertainty  can  be  detennined.  This  sensitivity  can 
serve  to  guide  further  exploration  in  the  reduction  of  model  uncertainty  through  additional 
information — specifically  detailing  the  need  for  posterior  model  likelihood  updating. 

While  model-form  uncertainty  quantifies  the  discrepancies  between  models,  parametric 
uncertainty  quantifies  the  difference  between  a  model  and  experimental  data.  Predictive 


64 


uncertainty  is  a  result  of  the  simplifying  assumptions  that  are  made  in  the  construction  of  a 
model,  such  as  an  inviscid  or  incompressible  flow  assumption.  Although  individually,  the 
models  are  deterministic,  the  presence  of  predictive  uncertainty  dictates  that  they  should  be 
instead  described  as  distributions,  to  account  for  the  known  errors  as  a  result  of  simplifying 
assumptions.  To  transform  the  deterministic  solution  to  each  model  into  a  distribution,  an 
assumption  is  made  that  each  model’s  estimation  contains  a  residual  that  is  identically, 
independently  and  normally  distributed  (IDD),  (Equation  (35)). 


•A,  =  y,  -  fk  (a  )  ~  N(0,  ak )  (35) 

In  Equation  (35)  yi  represents  an  experimental  result  at  design  variable  vector  xi  while 

fk  (xi )  represents  model  k’s  solution  at  the  same  design  variable  vector.  <jk  represents  the 
standard  deviation  of  the  nonnal  distribution,  as  calculated  by  Equations  (36): 


ZO'f  -fkixi)Y 

1=1 _ 

n 


(36) 


Now  that  the  residual  of  each  model  is  defined,  the  predictive  distribution  of  each  of  the 
models  can  be  constructed  by  adding  the  residual  to  the  detenninistic  model  prediction,  as  shown 
in  Equation  (37): 

p(y\Mk)=fk(x)  +  ek  (37) 
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Finally,  the  posterior  likelihood  of  each  model,  p(D  \  Mk),  can  be  calculated  by 


computing  the  joint  probability  of  the  experimental  data  as  shown  in  Equation  (38): 


p(D\Mk)  =  p(yl,...,y„  \Mk)  =  Y[p(yi\fk(xi),(jk)  (38) 

i=i 

The  posterior  likelihoods  for  each  model  calculated  by  Equation  (38)  can  then  be  used  to 
update  the  variable  responses  predicted  by  the  adjustment  factors  approaches.  This  technique 
allows  for  further  refinement  of  the  model  probabilities,  P{Mi),  if  deemed  necessary  by  the 

probabilistic  adjustment  factors  approach.  As  refinement  requires  the  addition  of  experimental 
data,  which  can  be  expensive  or  even  infeasible  at  an  early  design  stage,  it  is  crucial  to  utilize  the 
sensitivity  in  the  probabilistic  adjustment  factors  approach  to  estimate  the  merit  of  adding 
additional  data  into  the  calculation  of  the  model  likelihoods. 

To  demonstrate  the  application  of  model  uncertainty  to  an  aeroelastic  problem,  a  code 
was  written  that  solves  for  the  flutter  velocity  of  a  2  degree  of  freedom — pitching  and 
plunging — wing  subject  to  unsteady  aerodynamics  (Figure  31). 
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Figure  31:  2-DOF  Wing 


Including  the  presence  of  circulatory  flow,  the  lift  and  the  moment  about  the  shear  center 
and  the  moment  about  the  shear  center  can  be  calculated  by  Equations  (39)  and  (40)  respectively. 


Lsc  =  npb1\h-Ud-baci\+2npUbC(k') 


i  +  Ua  +  b\  —-a\d 


(39) 


Msc  =  npb 


(3ah  -Ub\  — a  j a-b1\^  +  cr  \ a 


+  2jipUb  \  a  +  -  \C(k) 


h  +  Ua  +  b\  —  -a\a 


(40) 


In  the  above  equations,  C(k )  represents  Theodorsen’s  Circulation  Function,  which 
controls  the  phasing  and  amplitude  of  the  lift  and  pitching  moments  with  respect  to  the  airfoil 
motion.  Theodorsen’s  Circulation  Function  is  a  complex  function  consisting  of  both  real  and 
imaginary  parts.  Multiple  surrogate  functions  exist  to  approximate  Theodorsen’s  Function  as  a 
function  of  k  over  the  range  of  k-values  experienced  by  the  system  (Equations.  (41  -  44)): 
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(41) 


C1(k)  =  1- 


0.165k  0.355k 

k  -0.0455/  k  -0.3/ 


0.01365 +  0.2808/k- 0.5k2 

CJk)  = - — 

0.01365  +  0.3455/A:  -k2 


C3(k)  = 


J\  J oyi 

J\  ~  J  oJ  J o  ~  J\  J 


C4(k)  = 


(1  +  10.61/A-X1  +  1.774/A-) 
(1  +  13.5 1/A:)(1  +  2.745/A:) 


Where  k  is  the  reduced  frequency  of  the  system  defined  by  Equation  (45)  where  to  is  the 
frequency  of  oscillation  of  the  airfoil  and  UK  is  the  free-stream  velocity 


co  *b 


The  real  and  imaginary  part  of  these  four  surrogate  models  vary  over  an  average 


operating  range  of  k’s,  as  shown  in  Figure  32. 
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Figure  32:  Real  (Left)  and  Imaginary(Right)  Components  of  C(k)  for  4  Models 


Once  a  surrogate  model  of  Theodorsen’s  Function  is  decided  upon,  the  flutter  velocity  for 
this  system  of  equations  can  then  be  solved  using  the  theory  of  unsteady  aerodynamics 
(Equations  (39  -  40))  and  the  V-g  Method.  However,  due  to  the  multiple  possible  surrogate 
models  to  use,  there  is  inherent  model  uncertainty  to  the  problem.  Executing  the  V-g  method  for 
the  sample  2  DOF  wing  produced  results  shown  in  Table  9. 


Table  9:  Flutter  Velocities  for  Four  Models 


Vf  (ft/s) 

Model  Probability 

Cl(k) 

101.393 

0.3 

C2(k) 

99.469 

0.15 

C3(k) 

97.968 

0.4 

C4(k) 

97.598 

0.15 

By  assigning  a  probability  to  each  of  the  four  models,  based  on  expert  opinion,  an 


adjustment  factors  approach  can  be  utilized  to  develop  a  distribution  of  the  flutter  velocity  for  the 
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2  DOF  system.  By  using  the  additive  and  multiplicative  adjustment  factors  approach  detailed  in 
the  sections  above,  a  normal  and  lognormal  distribution  of  the  flutter  velocity  can  be  developed. 
Table  10  shows  the  means  and  standard  deviations  of  the  figure  of  merit — in  this  case  the  flutter 
velocity  of  the  system — for  both  the  additive  and  multiplicative  adjustment  factors  approach. 
These  distributions  are  then  plotted  in  Figure  33,  showing  the  probability  density  functions  for 
the  flutter  velocity  of  the  wing  as  calculated  by  both  the  additive  and  multiplicative  adjustment 
factors  approach. 


Table  10:  Distribution  Parameters  for  Two  Methodologies 


Mean 

Standard  Dev. 

Additive  (normal) 

99.16515 

1.5635 

Multiplicative  (lognormal) 

99.15289 

1.55687 

Figure  33:  Normal  and  Lognormal  Plots  of  Vf 
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As  detailed  above,  the  primary  difference  between  multiplicative  and  additive  adjustment 


factors  approaches  is  the  assumption  of  the  final  distribution.  While  in  problems  where  the 
physical  modeling  is  well  understood  and  the  variance  between  models  is  relatively  small — such 
as  the  case  shown  above — the  two  approaches  might  result  in  similar  distributions,  in  cases 
where  the  initial  variance  between  models  is  large,  the  two  approaches  can  produced 
dramatically  different  result.  When  such  a  case  exists,  existing  knowledge  regarding  variable 
and  model  distributions  must  be  used  to  decide  upon  an  approach,  or  the  two  approaches  must  be 
analyzed  individually. 

The  next  step  is  to  begin  to  quantify  the  model  uncertainty  in  the  problem,  or  more 
specifically,  the  uncertainty  associated  with  the  probabilities  assigned  to  each  of  the  models.  By 
utilizing  the  probabilistic  adjustment  factors  approach  detailed  in  this  research,  the  model 
probabilities,  P(Mi )  ,  can  be  assigned  a  distribution  and  the  effect  of  their  inherent  uncertainty 

can  be  determined.  By  defining  each  of  the  four  model  probabilities  as  a  nonnal  distribution 
with  a  mean  of  their  deterministic  value,  shown  in  Table  10,  and  a  uniform  standard  distribution, 
the  effect  of  model  uncertainty  can  be  calculated.  With  a  standard  deviation  of  0.05  for  each 
model  probability  distribution,  the  means  and  standard  deviations  of  the  final  response  output  for 
flutter  velocity  were  calculated  and  can  be  seen  in  Table  11. 


Table  11:  Distribution  Parameters  for  Probabilistic  Adjustment  Factors  Approach 


Mean 

Standard  Dev. 

Additive 

99.1753 

1.4915 

Multiplicative 

99.1646 

1.4867 
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By  comparing  the  results  of  the  deterministic  and  the  probabilistic  adjustment  factors 
approach  (Table  10  and  Table  1 1),  it  can  be  seen  that  there  is  less  than  0.1%  percent  difference 
between  the  two  final  distribution  of  flutter  velocity  for  both  the  additive  and  multiplicative 
adjustment  factors  approach.  This  demonstrates  that  the  effects  of  model  uncertainty — the 
uncertainty  in  the  model  probabilities — are  minimal  for  this  problem,  and  that  additional 
quantification  of  predictive  uncertainty,  such  as  by  using  posterior  model  probability  updating, 
would  be  of  minimal  benefit  in  the  reduction  of  the  model  uncertainty.  This  is  valuable 
information  for  designers,  as  it  provides  information  regarding  the  relative  benefit  and  merit  of 
conducting  possibly  expensive  experiments  to  further  refine  and  reduce  model  uncertainty.  The 
results  of  this  study  would  show  to  the  designer  that  although  there  in  uncertainty  in  the  models 
being  used,  the  introduction  of  new  experimental  data  to  the  problem  would  not  further  reduce 
the  model  uncertainty. 
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6.0  OPTIMIZATION  SCHEME 


It  is  apparent  how  important  of  a  role  uncertainties  play  in  RLV  design  and  it  is  necessary 
to  reduce  the  amount  of  risk.  This  can  be  accomplished  by  introducing  an  optimization  scheme 
into  the  design  of  the  RLV  by  minimizing  risk.  In  most  cases  design  optimization  is  thought  to 
be  based  on  a  detenninistic  problem.  Typically  deterministic  design  optimization  fonnulates  the 
problem  with  an  objective  function  that  is  to  be  minimized  or  maximized  bounded  by  constraints 
which  can  be  seen  in  Equations  (46  -  48).  The  results  to  this  optimization  will  consist  of  a  single 
design  point  that  satisfies  the  constraint  conditions. 

Minimize: 


/(*) 


(46) 


Subject  to: 


giOO  ^  0 


(47) 


hj  (x)  =  0 


(48) 


In  this  case,  x  is  the  vector  of  design  variables  and  hj(x)  and  gi(x)  are  equality  constraints 
and  inequality  constraints  respectively. 
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Reliability-Based  Design  Optimization  (RBDO)  is  a  valuable  evolution  in  deterministic 
design  optimization.  Instead  of  optimizing  a  problem  based  on  a  deterministic  problem  the 
objective  function  is  based  on  the  probability  of  failure  found  using  uncertainty  quantification 
methods.  Reliability-based  design  optimization  accounts  for  these  statistical  distributions  in  the 
analysis  through  stochastic  finite  elements  using  FORM,  SORM,  evidence  theory  or  many  other 
reliability-based  methods.  RBDO  fonnulates  the  problem  with  an  objective  function  that  is  to  be 
minimized  or  maximized  bounded  by  constraints  which  can  be  seen  in  Equations  (49  -  5 1). 

Minimize: 


(49) 


Subject  to: 


giOO  ^  0 


(50) 


hj  (x)  =  0 


(51) 


In  this  case,  x  is  the  vector  of  design  variables  that  can  contain  uncertainties  or  can 
influence  uncertainties  depending  on  how  the  problem  is  established.  h,(x)  and  gi(x)  are  the  same 
as  in  the  detenninistic  design  optimization  problem.  The  objective  function  is  now  to  reduce  the 
probability  of  failure  resulting  in  a  highly  reliable  design. 
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Now  that  a  basic  understanding  of  reliability-based  design  optimization  is  understood  an 
outline  of  the  proposed  design  optimization  of  the  RLV  wing  based  on  aeroelastic  flutter 
simulation  is  explained.  Due  to  the  wide  presence  of  uncertainty,  as  well  as  the  different  types  of 
uncertainty  present,  accounting  for  uncertainty  in  the  problem  becomes  necessary  to  perform  a 
complete  design.  Figure  34  illustrates  a  flow  chart  of  the  propagation  of  uncertainties 
incorporated  in  an  RBDO  algorithm.  In  the  RBDO  section  three  methods  will  be  explored; 
gradient-based  RBDO,  practical  swann  optimization,  and  a  proposed  method  of  adaptive  particle 
swarm  optimization. 
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Figure  34:  Vehicle  Design  Environment  Framework 
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Design  Optimization  Design 


6.1  Reliability-Based  Design  Optimization  Incorporating  Evidence  Theory 

Reliability-Based  Design  Optimization  (RBDO)  techniques  are  developed  to  address  the 
analytical  guarantee  of  the  performance  of  a  structural  system.  In  this  section,  uncertainty 
quantification  using  evidence  theory  demonstrated  in  the  previous  sections  is  implemented  into 
design  optimization.  To  address  the  discontinuity  of  the  measurements  (BEL  and  PL),  a 
supplementary  measurement,  plausibility  decision,  described  in  section  5.3,  will  be  used.  The 
flow  chart  of  gradient-based  design  optimization  is  demonstrated  in  Figure  35. 


Figure  35:  Gradient  RBDO  Flow  Chart 
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As  stated  in  section  5.3.1  the  four-bar  truss  structure  (Figure  36)  will  be  revisited  to 
investigate  an  RBDO  analysis.  To  complete  the  design  optimization  portion  of  the  analysis 
gradient  information  is  needed  to  approximate  the  constraint  function  because  no  closed  form 
equation  can  be  derived  when  evidence  theory  is  implemented.  A  surrogate  model  of  the 
constraint  was  developed  using  TANA  approximation  (explained  in  section  5.2.1.)  This 
approximation  is  only  valid  in  short  intervals  the  process  outline  in  Figure  35  and  requires 
iterations  to  get  an  accurate  representation  of  the  function. 

The  problem  fonnulation  for  this  example  is  as  follows: 

Minimize: 


/  =  3x1  +  V3x2 


(52) 
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Subject  to: 


9 1 


-^>0 

x2 


(53) 


92=*1~  5-73  ^  0 


(54) 


93  =  x2  -  7.17  >  0  (55) 

drf  —  Pldec  9 1  success  >  0  (56) 

Where  drf  is  the  desired  reliability  factor  for  the  constraint  gi  drf  is  chosen  by  a 
designer  based  on  how  reliably  they  want  the  structure.  Table  12  contains  the  results  from  the 
gradient-based  RBDO  using  Benanazer’s  plausibility  decision  approximation.  As  seen  in  the 
table  there  are  a  range  of  optimizations,  in  which  the  drf  was  altered.  The  results  show  that 
since  this  problem  was  constrained  with  PL  dec  rather  than  BEL  or  PL  that  the  desired 
reliability  factor  was  achieved  in  each  analysis.  Notice  as  the  drf  is  increased  (meaning 
reliability  of  the  structure  is  increasing)  the  higher  the  objective  function  value  is.  This 
demonstrates  a  more  reliable  structure  cost  more  materials  then  a  less  reliable  one. 
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Table  12:  Four-Bar  Truss  Reliability  Design  Optimization  Results 


f 

xl 

x2 

drf 

PL_dec 

BEL 

PL 

43.52313 

9.19752 

9.19752 

0 

0 

0 

0.015 

43.99996 

9.335168 

9.234402 

0.01 

0.01 

0 

0.0255 

44.90086 

9.423517 

9.601514 

0.1 

0.1 

0.015 

0.253 

45.26616 

9.627671 

9.458821 

0.2 

0.2 

0.055 

0.6095 

45.44449 

9.541425 

9.711154 

0.3 

0.3 

0.175 

0.6095 

45.62129 

9.468467 

9.939728 

0.4 

0.4 

0.175 

0.6195 

45.78633 

9.580633 

9.839363 

0.5 

0.5 

0.175 

0.91 

46.03484 

9.682678 

9.807335 

0.6 

0.6 

0.53 

0.917 

46.35316 

9.780256 

9.822107 

0.7 

0.7 

0.54 

0.917 

46.67998 

9.785824 

10.00115 

0.8 

0.8 

0.54 

1 

47.18891 

9.962875 

9.988322 

0.9 

0.9 

0.9 

1 

115.6673 

23.78244 

25.58815 

0.95 

0.95 

0.9 

1 

152.1028 

27.18551 

40.57945 

0.99 

0.99 

0.9 

1 

159.0244 

26.79981 

45.39415 

1 

1 

1 

1 

6.2  Particle  Swarm  Optimization 


Vehicle  design  is  an  inherently  non-linear  multi  physical  problem  with  the  presence  of 
continuous,  mixed,  and  integer  variables.  To  analyze  the  relative  merit  of  a  particular  design, 
some  degree  of  simulation  required,  and  the  time  required  for  this  simulation  varies  with  the 
fidelity  of  the  analysis  required.  In  addition,  the  detennination  of  design  variable  gradients  is  not 
analytically  possible  when  simulation  beyond  basic  static  finite  elements  is  employed. 

Traditional  gradient-based  optimization  methods  rely  on  the  input  of  design  variable  gradients, 
which  would  need  to  be  detennined  using  a  finite-differencing  scheme  or  analytical  gradients — 
which  are  unavailable  or  difficult  to  obtain  for  many  complex  simulation  analyses.  Due  to  the 
relatively  long  simulation  time  required  for  an  analysis,  utilizing  a  finite  differencing  scheme  can 
quickly  become  computationally  restrictive.  In  addition,  gradient-based  methods  perfonn  very 
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poorly  in  highly  non-linear  problems,  finding  only  local  optima  simply  as  a  function  of  the  initial 
design  point. 

Traditional  heuristic  methods,  such  as  particle  swann  optimization  (PSO)  help  to 
alleviate  some  of  these  design  problems  [45].  PSO  (Equations  (57,  58))  is  a  heuristic  method 
that  will  converge  upon  a  global  optimum  and  does  not  require  gradient  information  at  particular 
points. 


(57) 


xlk+ 1  =  xik  +vik+1At 


(58) 


However,  as  with  most  heuristic  algorithms,  PSO  requires  numerous  function 
evaluations — often  within  a  very  small  region  of  the  design  space — to  converge  upon  an  optimal 
point.  In  addition,  due  to  the  large  number  of  user-controlled  parameters  in  the  optimization,  the 
algorithm  can  perfonn  much  differently  depending  on  the  values  of  those  inputs,  which  are  based 
on  the  user’s  knowledge  of  the  system  at  hand.  Due  to  the  lengthy  simulation  time  required  in 
vehicle  design  evaluations;  this  once  again  becomes  computationally  restrictive.  When 
uncertainties  within  the  design  problem  are  introduced,  the  computational  cost  becomes  even 
more  restrictive,  and  traditional  optimization  methods  become  costly  to  the  point  of  irrelevance. 

In  order  to  optimize  the  design  of  a  proposed  vehicle  or  vehicle  component,  these 
uncertainties  listed  above  must  be  accounted  for  and  an  efficient  algorithm  must  be  able  to 


81 


converge  upon  an  optimum  in  an  acceptable  amount  of  time.  With  the  introduction  of 
uncertainties  into  the  problem,  the  number  of  function  evaluations  required  can  increase 
dramatically.  Ko  et  al  [46]  has  previously  used  particle  swann  to  solve  a  trajectory  optimization. 
Dimou  et  al  also  explored  reliability-based  optimal  design  of  truss  structures  using  particle 
swarm  optimization  [47].  Thus,  a  method  must  be  implemented  for  the  optimization  routine  that 
minimizes  the  number  of  function  evaluations  required  while  still  converging  upon  an  optimal 
solution. 

The  problem  executed  with  this  study  is  a  closed  form  problem  to  stand  as  an  initial 
surrogate  for  an  actual  vehicle  design  problem.  This  was  done  to  limit  the  computational  time 
required  for  an  initial  study,  as  the  number  of  simulations  required  for  each  optimization  can 
reach  the  thousands  and  a  complete  initial  design  space  exploration  was  executed — meaning  that 
nearly  one  million  simulations  would  have  been  required  for  this  initial  study.  The  optimization 
problem  executed  for  this  was  the  “Egg  Crate  Function”  (Equation  (59));  a  highly  non-linear 
problem  with  multiple  optima  in  a  small  region  that  is  often  used  in  the  literature  to  explore  the 
relative  merit  of  optimization  routines  in  highly  non-linear  design  environments. 

Minimize: 


f{x)  =  xt2  +  x22  +  25(sin(x1)2  +  sin(x2)2)  (59) 
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Subject  to: 


— 2n  <  x1  <  2tt 


(60) 


— 2tt  <  x2  <  2tt 


(61) 


To  initiate  the  design  space,  Latin  Hypercube  Sampling  was  used  to  span  the  entire 
design  space  with  initial  sampling  points  and  ensure  faster  convergence.  Figure  37  below  shows 
a  plot  of  the  function  and  a  sample  iteration  history  for  one  of  the  optimizations  of  the  “Egg 
Crate  Function”. 
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Figure  37:  “Egg  Crate  Function”  and  Sample  Iteration  History 
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For  this  closed-form  optimization  problem,  a  design  space  study  was  conducted  to 
detennine  the  number  of  function  evaluations — which  directly  correlate  to  the  number  of 
simulations  required  to  analyze  a  vehicle’s  performance — to  converge  upon  an  optimal  solution. 
This  included  a  study  on  the  effect  of  the  different  parameters  in  particle  swarm  optimization  to 
see  how  they  affect  the  efficiency  of  the  algorithm.  For  this  study,  ten  different  optimizations 
were  executed  at  a  baseline  parameter  set  to  establish  an  initial  performance  value  Table  13. 

Then,  individually,  each  controllable  parameter  was  changed  though  the  full  spectrum  of  values 
to  determine  which  parameters  had  significant  effects  on  the  efficiency  of  the  algorithm.  Plots  of 
the  number  of  simulations  required  for  convergence  versus  the  change  in  parameter  values  can 
be  seen  in  Figure  38,  Figure  39,  Figure  40,  and  Figure  41.  In  addition,  for  the  purpose  of  this 
initial  study,  only  first-order  variable  interactions  were  explored. 


Table  13:  “Egg  Crate  Function”  Optimization  results  at  standard  parameters 


Num_part 

Num_neigh 

fw 

w 

cl 

c2 

yo 

No.  Fun.  Calls 

20 

5 

0.975 

0.9 

1.5 

1.5 

4.221  E-08 

940 

20 

5 

0.975 

0.9 

1.5 

1.5 

4.158E-05 

960 

20 

5 

0.975 

0.9 

1.5 

1.5 

8.126E-05 

900 

20 

5 

0.975 

0.9 

1.5 

1.5 

7.949E-06 

1000 

20 

5 

0.975 

0.9 

1.5 

1.5 

6.730E-04 

940 

20 

5 

0.975 

0.9 

1.5 

1.5 

3.454E-04 

1060 

20 

5 

0.975 

0.9 

1.5 

1.5 

3.340E-05 

1040 

20 

5 

0.975 

0.9 

1.5 

1.5 

5.925E-04 

980 

2.219E-04 

977.5 
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Neighboring  Particles  vs.  Simulations 
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Figure  38:  Number  of  Neighbors  vs.  Function  Evaluations 


C-Values  vs.  Simulations 


Cl  &  C2  Values 


Figure  39:  Individuality  and  Socialability  Factors  vs.  Function  Evaluations 
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Number  of  Particles  vs.  Simulations 


Number  of  Particles 


Figure  40:  Number  of  Particles  vs.  Function  Evaluations 


Convergence  Value  vs.  Simulations 
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Figure  41:  Convergence  Value  vs.  Function  Evaluations 
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It  can  be  seen  from  these  results  that  with  traditional  particle  swann  optimization,  some 
of  the  input  parameters  into  the  algorithm  have  dramatic  effects  upon  the  number  of  functions 
required  for  convergence.  In  addition  to  the  data,  the  quality  of  the  solutions  that  the  algorithm 
converged  to  with  different  parameter  values  was  recorded,  but  it  was  found  that  the  quality  of 
the  solution  was  only  a  function  of  the  convergence  criteria,  and  not  of  any  additional  parameters 
in  the  problem. 

After  identifying  significant  input  parameters  in  the  conventional  particle  swann 
algorithm,  a  modification  was  made  to  the  algorithm  in  attempts  to  reduce  the  number  of 
simulations  required  for  convergence  while  still  maintaining  the  desired  level  of  accuracy.  This 
method,  called  Adaptive  Particle  Swarm  Optimization  (A-PSO),  is  detailed  in  section  6.3. 

6.3  Adaptive  Particle  Swarm  Optimization  (A-PSO) 

As  shown  in  the  previous  section,  the  computational  cost  for  particle  swarm  optimization 
can  be  great.  For  a  simple  two-dimensional  problem,  it  was  shown  that  up  to  1,000  simulations 
could  be  required  for  convergence,  and  this  is  before  any  degree  of  uncertainty  is  introduced  into 
the  problem — which  as  shown  in  previous  sections  could  increase  the  number  of  simulations  by 
a  factor  of  20  for  a  simple  two  dimensional  problem.  In  addition,  as  greater  dimensionality  is 
added  to  the  problem,  these  numbers  scale  exponentially;  resulting  in  uncertain  designs  with  only 
5-6  design  variables  requiring  upwards  of  1,000,000  simulations  for  convergence  using 
traditional  particle  swarm  optimization.  With  detailed  simulations,  this  number  of  required 
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simulations  becomes  prohibitively  costly,  and  an  alternative  method  must  be  found  to  reduce  this 
number  of  simulations  to  a  more  manageable  number. 

While  particle  swarm  optimization  requires  a  large  number  of  function  evaluations  to 
converge  upon  an  optimal  solution,  it  can  be  seen  that  a  large  number  of  these  function 
evaluations  occur  within  a  small  region  of  the  design  space.  By  implementing  surrogate  models 
in  these  small  regions  of  the  design  space,  the  number  of  actual  function  evaluations  can  be 
dramatically  reduced.  However,  when  using  surrogate  modeling,  there  is  an  added  error  due  to 
the  inaccuracy  of  the  model  being  implemented.  In  addition,  if  these  models  are  considered 
valid  over  large  regions  of  the  design  space,  the  error  is  increased.  Thus,  it  can  be  seen  that  the 
small  region  of  the  design  space  that  a  surrogate  model  is  considered  valid,  the  more  accurate  it 
generally  is.  However,  if  the  entire  design  space  is  simply  partitioned  into  multiple  small  regions 
to  be  modeled  with  surrogate  models,  the  number  of  simulations  required  to  fonn  these  models 
would  once  again  be  large  and  unwieldy.  Thus,  what  was  done  was  to  initiate  traditional  particle 
swarm  optimization.  The  design  space  is  then  divided  into  several  large  regions  to  be  modeled 
using  a  second  order  response  surface  model.  Equation  (62) 

Yj=Po+  HU  frXji  +  Y.U  PiXji 2  +  Hi<  Hnv  Pip  Xtj  XJp  +  sp  (62) 

For  the  next  iteration  in  the  optimization,  the  function  evaluations  can  then  be 
approximated  by  the  response  surface  models,  which  are  considered  valid  within  their  region. 
Then,  in  order  to  add  additional  data  into  the  models,  a  finite  number  of  points  (an  additional 
parameter  added  to  the  traditional  PSO  algorithm)  are  randomly  selected  for  full  evaluation. 
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These  points  are  fully  evaluated — in  the  closed  fonn  problem  this  is  rather  simple,  but  in  a 
simulation-based  problem,  this  would  represent  additional  simulations — and  this  new  data  is  then 
used  in  conjunction  with  the  existing  points  to  update  the  surrogate  models.  When  a  region  of 
the  design  space  defined  by  a  model  becomes  over-populated  with  function  evaluations,  the 
space  is  then  partitioned  into  smaller  regions,  reducing  the  bounds  on  the  models  within  that 
region.  This  can  be  best  illustrated  in  Figure  42. 


Figure  42:  Model  Design  Space  Reduction  in  A-PSO 


This  design  space  reduction  is  designed  so  that  as  the  optimization  algorithm  converges 
upon  a  solution,  the  surrogate  models  become  smaller  and  more  numerous  in  the  area  of  multiple 
particles,  which  corresponds  to  the  location  of  an  optimal  solution.  Thus,  in  this  method, 
previous  data  regarding  function  evaluations  is  not  wasted,  but  used  to  form  surrogate  models  to 
aid  in  the  convergence  of  the  algorithm. 
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7.0  FAST  WING  RESULTS 


7.1  Aleatory  Uncertainty  Investigation 

To  further  investigate  FORM  and  SORM  for  RLV  design  an  aeroelastic  model  was 
analyzed  with  structural  aleatory  uncertainties.  Using  the  model  described  in  the  chapter  4  an 
aeroelastic  flutter  uncertainty  quantification  analysis  was  investigated  (Figure  43). 


Figure  43:  NASTRAN  FEA  Wing  Model  for  Aleatory  Risk  Quantification  Analysis 


The  two  uncertain  variables  that  were  included  in  this  analysis  were  selected  to  be  the 
skin  thickness  and  cross-sectional  area  of  the  wing.  These  variables  are  uncertain  due  to  complex 
manufacturing  process  of  composite  materials.  The  parameters  used  for  the  atmospheric 
conditions  in  this  analysis  represent  the  vehicle  entering  the  atmosphere  during  the  trajectory  at 
subsonic  levels.  The  limit  state  allowed  the  flutter  dynamic  pressure  to  reach  2000  psf  before  it 
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was  considered  as  a  failure.  The  uncertain  variables  were  represented  by  nonnal  distributions. 
The  means  and  standard  deviations  of  the  skin  and  spars  cross-sectional  area  were  p=.04  c=.005 
and  p=.2  c=.02  inches  respectively.  The  results  for  the  aleatory  uncertainty  quantification  of  the 
aeroelastic  model  are  located  in  Table  14. 


Table  14:  Aleatory  Uncertainty  Quantification  Results  Including  Aeroelastic  Analysis 


Total  Number  of 
Function  Calls 

Total  Computational 
Time 

Pr 

%  Difference 
from  MC 

FORM  with  adaptive 
approximations 

18 

22  seconds 

0.05075 

.132 

Breitung's  Fonnulation 

23 

29  seconds 

0.05082 

.006 

Tvedt's  Formulation 

23 

29  seconds 

0.05084 

.045 

Koyluoglu's  Formulation 

23 

29  seconds 

0.05066 

.309 

Monte  Carlo 

150,000 

33.33  hours 

0.05082 

0 

7.2  Epistemic  Uncertainty  Investigation 

This  study  focuses  on  a  reliability  analysis  of  an  RLV  flutter  speed.  In  past  aerospace 
designs  the  vehicles  were  modeled  without  considering  uncertainty  which  may  have  led  to  the 
Space  Shuttle  Challenger  and  the  more  recent  Space  Shuttle  Columbia  disasters.  To  quantify  the 
uncertainties  a  more  robust  design  can  be  implemented  to  the  space  program  to  avoid 
catastrophic  events. 
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Epistemic  uncertainty  is  a  type  of  uncertainty  where  little  infonnation  is  known  regarding 
the  variable,  and  it  would  be  inaccurate  to  assume  some  sort  of  distribution.  Only  epistemic 
uncertainties  will  be  considered  in  this  investigation.  This  problem  will  be  formulated  using  non 
probabilistic  methods  since  evidence  theory  is  being  incorporated.  The  uncertain  variables  are 
selected  in  this  problem  based  on  the  little  infonnation  known  about  the  variable  and  the 
importance  they  play  in  the  vehicle  design.  For  example  the  composite  skin  on  the  RLV  where 
two  samples  of  the  same  composite  can  have  completely  different  properties.  This  is  related  to 
the  orientation  of  the  layers.  There  is  not  enough  infonnation  to  construct  a  pdf  on  the  orientation 
of  each  thus  the  variable  is  epistemic.  Since  RLV  design  consist  of  so  many  types  of  variables  a 
limit  state  must  be  declared  in  order  to  complete  an  uncertainty  analysis. 

Flutter  will  be  used  as  the  limit  state  of  this  problem.  In  complex  structures  where  both 
the  aerodynamics  and  the  mechanical  properties  of  the  structure  are  not  fully  understood  the 
aeroelastic  flutter  phenomenon  described  in  the  previous  section  can  be  quantified  with 
uncertainties.  As  seen  in  Figure  44  a  demonstration  of  an  uncertainty  bound  is  added  to  the 
flutter  assessment.  This  bound  shows  that  if  uncertainty  is  incorporated  to  an  analysis  it  can  cross 
the  flutter  margin  causing  catastrophic  failure. 
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5000 


Flutter  Dynamic  Pressure  Symmetric  Boundary  Condition 


Figure  44:  Aeroelastic  Analysis  with  Uncertainty  Bound 


The  variables  that  were  chosen  for  this  problem  are  based  on  two  parts  of  the  aeroelastic 
analysis  to  give  diversity  to  the  problem  to  demonstrate  the  coupling  of  uncertainties  between  the 
aerodynamic  and  structural  model.  The  uncertainties  come  from  the  composites  in  the  structural 
model  and  the  atmospheric  conditions  in  the  aerodynamic  model.  As  seen  in  Figure  45  on  the  left 
is  the  makeup  of  the  composite  material  for  the  skin  of  the  RLV  and  the  right  is  the  equation  to 
calculate  the  air  density.  Composites  have  many  uncertainties  associated  with  them  because  they 
are  hard  to  reproduce  consistently.  The  uncertainties  in  the  composite  for  this  particular  problem 
come  from  the  thickness  which  can  vary  of  the  composite  as  well  as  the  orientation  of  the  third 
and  fourth  layer.  These  are  considered  epistemic  uncertainties  because  there  is  not  enough 
information  on  this  type  of  composite  to  construct  a  pdf.  The  second  group  of  uncertainties 
comes  from  the  air  density.  The  three  uncertain  variables  associated  with  air  density  are  the  air 
pressure,  gas  constant,  and  temperature.  The  gas  constant  variable  was  selected  to  simulate 
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moisture  in  the  atmosphere.  The  temperature  and  pressure  were  chosen  because  air  pressure  has 
a  large  dependence  on  them.  The  atmospheric  conditions  are  thought  to  be  epistemic  uncertainty 
because  there  is  not  enough  infonnation  to  accurately  construct  a  pdf  without  introducing  more 
uncertainty  into  the  problem. 


Layer  1 
Layer  2 
Layer  3 
Layer  4 


Figure  45:  Uncertainties  for  six  variables 


Figure  46  represents  the  three  uncertain  variables  BBA’s  for  the  composite  uncertainties 
and  Figure  47  signifies  the  three  uncertain  variables  BBA’s  in  the  atmospheric  uncertainties.  In 
this  problem  three  experts  were  asked  their  opinion  on  each  variable,  at  which  interval  they 
thought  the  most  likely  occurrence  of  uncertainty  would  occur.  Each  interval  was  then  weighted 
on  the  expert’s  opinion.  Figure  46  and  Figure  47  illustrates  each  expert’s  opinion  in  which  for 
the  variables  layer  thickness  and  temperature  all  three  opinions  were  used.  In  the  other  four 
variables  one  expert  did  not  express  enough  knowledge  in  that  particular  subject  and  the  variable 
was  omitted. 
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Composite  Uncertainties 
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Figure  46:  Composite  Uncertainties  Basic  Belief  Assignments 
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Figure  47:  Aerodynamic  Uncertainties  Basic  Belief  Assignments 
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As  mentioned  in  the  two  variable  demonstration  case,  the  limit  state  for  the  six  variable 
problem  is  also  2000  psf.  The  results  for  the  six  variable  problem  are  depicted  in  Figure  48.  The 
reliability  of  the  RLV  falls  between  the  bounds  BEL  and  PL.  This  is  expected  since  the  RLV  is 
initially  designed  to  have  a  high  relibility.  The  bounds  can  be  improved  to  find  a  more  accurate 
reliability  and  are  explored  in  the  next  study. 


BEL(A)  =  V  m(C) 


PL(A)  =  I  m(C) 


CcA  CnA^O 

BEL  =  0.925  PL  =1.0 

Figure  48:  Six  Variable  RLV  Problem  Results 


_ Simulations _ 3456 

CPU  time  per  simulation  30  sec 

Total  Time  28.8  hours 

Figure  49:  Simulation  Time  Six  Variable  RLV  Problem 


Figure  49  illustrates  the  efficiency  of  evidence  theory.  Although  the  method  makes  no 
assumptions  and  gives  an  accurate  bound  of  the  limit  state  the  computational  time  is  expensive 
The  more  variables  with  uncertainty  the  more  complex  the  problem  becomes.  As  seen  in  the 
figure  it  took  28.8  hours  to  execute  the  analysis  when  each  of  the  simulations  took  30  sec  to 
complete.  This  is  why  a  model  must  have  a  short  execution  time  otherwise  the  analysis  would 
take  years  to  complete.  One  possible  solution  would  be  to  make  a  surrogate  model  of  the 
aeroelastic  flutter  analysis  such  as  a  response  surface  which  would  look  act  like  a  closed  form 
solution,  resulting  in  a  lower  computational  time.  The  problem  with  a  surrogate  model  is  it 
would  introduce  uncertainty  to  the  problem  which  is  trying  to  be  avoided  when  using  evidence 
theory. 
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The  pay  off  of  using  evidence  theory  compared  to  a  probabilistic  method  is  no  new 
uncertainty  was  added  to  the  problem  by  assuming  a  pdf  for  the  uncertain  variables.  If  the  wrong 
pdf  was  assumed  and  a  probabilistic  method  such  as  first  order  reliability  method  or  second  order 
reliability  method  was  used  the  reliability  of  the  analysis  could  have  deviated  from  the  actual 
reliability  causing  misleading  infonnation. 

There  are  two  ways  to  improve  the  results  of  this  type  of  analysis  one  being  improving 
the  BBA  structures  and  the  second  by  introducing  the  concept  of  Plausibility  decision  (PL  dec). 
If  the  plausibility  belief  bound  is  too  large  one  solution  is  to  further  discretize  the  problem.  This 
can  be  accomplished  by  improving  the  BBA  structures.  To  improve  the  BBA’s  one  can  ask  more 
experts,  complete  more  experiments,  or  assume  the  expert  opinions  are  evenly  distributed  and 
split  them  up  into  more  pieces.  Improving  the  BBA’s  structure  will  shrink  the  bound  of  the  belief 
and  plausibility  giving  a  more  understood  problem. 

One  concern  for  an  analysis  with  flutter  is  what  if  flutter  does  not  occur  during  one  of  the 
iterations.  Evidence  theory  does  not  use  iterations  so  one  simulation  does  not  depend  on  the  next. 
This  works  out  very  nicely  when  a  simulation  does  not  flutter.  It  can  be  considered  as  a 
successful  case  and  the  analysis  could  continue  without  penalty.  For  the  current  demonstration 
flutter  occurred  at  every  simulation 
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7.3  Reliability-Based  Design  Optimization  Epistemic  Uncertainty  Investigation 

Now  that  evidence  theory  has  been  demonstrated  in  the  previous  section  it  is  apparent  the 
probability  of  success  can  be  improved.  In  this  study  an  RLV  wing  structure  will  be  optimized  to 
reduce  the  weight  while  increasing  the  probability  of  success  based  on  a  flutter  limit  state. 

The  first  step  in  the  problem  is  to  identify  uncertainties.  As  mentioned  before  composite 
materials  have  uncertain  properties  and  are  hard  to  manufacture  consistently.  In  the  previous 
section  the  uncertainties  came  from  the  orientation  of  the  layers  of  the  composites  that  make  up 
the  skins  of  the  wings.  In  this  study  only  the  composite  thickness  on  the  wings  skin  will  be 
examined.  Figure  50  contains  images  of  the  three  sets  of  skin  thicknesses  that  are  being 
examined.  The  blue  section  will  be  referred  to  as  the  wing  tip  skin,  the  red  portion  will  be 
recognized  as  the  top  skin,  and  the  green  segment  will  be  identified  as  the  bottom  skin. 


98 


Figure  50:  Variability  Design  of  Skin  Thickness  Variables 

This  study  also  includes  atmospheric  density  as  an  uncertainty.  Like  before  the 
atmospheric  density  is  an  epistemic  variable  because  there  is  not  enough  information  to 
detennine  a  well  defined  pdf.  Instead  of  splitting  the  air  density  into  an  equation  and  three 
separate  uncertain  demonstrated  in  the  previous  section  it  will  consist  of  only  one  uncertainty. 
This  was  done  to  reduce  the  computational  time  because  now  the  problem  is  an  iterative  process 
and  each  additional  uncertain  variable  increases  computational  time  dramatically  when 
employing  evidence  theory.  The  BBAs  for  the  four  uncertain  variables  can  be  found  in  Figure 
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51.  Like  before  in  the  RBDO  using  evidence  theory  example  the  three  skin  thicknesses  are 
distributed  in  percentages  because  not  only  will  the  skin  thicknesses  be  the  uncertainties  in  the, 
they  will  also  be  the  design  variables  in  which  the  problem  is  optimized. 


Wing  Tip  Skin  Thickness 

90%  92%  94%  96%  98%  100%  102%  104% 

< - X - X - > 

.16  .39  .45 

Wing  Top  Skin  Thickness 

84%  86%  88%  90%  92%  94%  96%  98%  100%  102% 

< - X - X - > 

.08  .75  .16 

Wing  Bottom  Skin  Thickness 

92%  94%  96%  98%  100%  102%  104%  106% 

< - X - X - > 

.12  .82  .06 

Air  Density 

1.00E'7  1.05E'7  1.10E'7  1.15E7  1.20E'7  1.25E'7  1.30E7  1.35E'7  1.40E'7 

< - X - X - > 

.17  .63  .20 

Figure  51:  BBAs  for  the  Wing  structure  used  in  the  RBDO  study 


The  BBAs  were  developed  based  on  three  expert  opinions.  Like  in  the  evidence  theory 
study  of  the  RLV  wing  the  limit  state  of  the  uncertainty  quantification  algorithm  is  flutter.  While 
completing  the  evidence  theory  analyses  during  the  RBDO,  if  a  flutter  dynamic  pressure  exceeds 
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2000  psf,  it  is  considered  a  failure.  The  goal  of  the  optimization  is  to  minimize  the  weight 
(objective  function.  Equation  (63))  of  the  wing  while  keeping  a  plausibility  decision  above  .98. 
Leading  to  the  problem  being  defined  as: 


Minimize: 


w  =  l,UPiVi 


(63) 


Subject  to: 


gt  =  flutter  limit  state(2000psf)  >  0  (64) 


drf  —  Pldec  9 1  success  >  0  (65) 

Where  n  is  the  number  of  elements  that  consist  of  the  skin  section.  pt  and  Vi  are  the 
mass  density  and  volume,  respectively,  of  the  ith  structural  element  participating  in  the 
design,  drf  is  the  desired  reliability  factor  which  in  this  case  is  .98. 


Table  15  shows  the  initial  and  optimized  results  for  the  RLV  wing  optimization.  The 
optimization  resulted  in  an  increase  in  PL  dec  and  Belief  while  reducing  the  weight  of  the  wing. 
The  optimization  converged  at  15  iterations.  Figure  52  shows  the  computational  time  invested  in 
this  RBDO. 


Table  15:  RBDO  Wing  Thickness  Results 


Wing-Tip 
skin  (in) 

Wing-top 
skin  (in) 

Wing-Bottom 
skin  (in) 

Weight 

(lb) 

PL_dec 

PL 

BEL 

Initial 

0.0576 

0.144 

0.0536 

769.78 

0.9175 

0.99 

0.786 

Optimized 

0.0488 

0.132 

0.0632 

766.5072 

0.98 

0.99 

0.8242 
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15 


_ Iterations 

CPU  time  per  Iteration  8  hours 
Total  time  120  hours 

Figure  52:  Simulation  Time  RBDO  Analysis 
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8.0  DISCUSSION  AND  CONCLUSIONS 


The  technical  effort  conducted  in  this  task  developed  methods  relevant  to  the  reliability- 
based  structural  design  of  an  RLV.  The  first  step  completed  in  the  analysis  of  the  RLV  wing  was 
a  trajectory  optimization  for  a  rocket-back  mission,  where  it  was  discovered  that  critical  points  in 
the  trajectory  needed  further  investigation.  From  the  mission’s  trajectory,  critical  points  were 
selected  based  on  maneuverability  and  rapid  changes  in  atmospheric  conditions. 

To  evaluate  these  critical  points,  a  finite  element  model  was  needed.  The  wing  model  was 
based  on  representative  dynamic  characteristics  of  the  RLV  configuration.  Since  no  physical 
model  of  the  RLV  wing  is  available,  the  created  wing  model  was  validated  based  on  the  original 
model  using  two  methods  of  characteristics  validation.  The  two  methods  used  were  the  modal 
assurance  criterion  for  mode  shape  comparisons  and  frequency  comparison  analysis.  From  these 
analyses  we  concluded  that  the  two  models  compared  are  in  sufficient  agreement  that  the  created 
model  can  be  used  as  a  surrogate.  Along  with  the  structural  model  an  aeroelastic  model  was 
constructed  in  which  a  series  of  flutter  analyses  were  conducted.  These  analyses  were  based  on  a 
flight  envelope  during  the  launch  phase  of  the  trajectory.  This  regime  was  selected  because  the 
RLV  experienced  the  most  dynamic  pressure  during  this  phase.  From  this  study  it  was 
detennined  that  the  critical  Mach  number  relating  to  the  flutter  speed  was  at  Mach  1.1. 

Once  the  model  was  validated  and  finalized,  uncertainty  quantification  methods  were 
explored.  Three  sources  of  uncertainty  were  explored  in  this  research:  aleatory,  epistemic,  and 
model-form.  It  was  discovered  that  when  quantifying  aleatory  uncertainties,  the  second-order 
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reliability  method  produces  better  results  than  the  first-order  reliability  method  at  the  cost  of 
more  computational  time.  Epistemic  uncertainty  was  analyzed  using  evidence  theory  where 
uncertain  variables  with  limited  infonnation  were  quantified.  This  resulted  in  bounds  of  actual 
response.  A  new  metric  known  as  plausibility  decision  was  introduced  to  estimate  the  reliability 
within  the  bound  and  to  obtain  gradient  information.  Three  methods  of  plausibility  decision  were 
investigated,  and  Benanzer’s  approximation  was  determined  to  be  the  most  fitting  for  an 
evidence  theory  analysis  with  more  than  two  variables.  This  effort  also  produced  a  new  method 
to  quantify  model-form  uncertainty.  This  type  of  uncertainty  is  present  when  more  than  one  type 
of  model  is  available  for  an  analysis,  such  as  finite  element  models  with  different  fidelities  or 
using  different  computational  methods. 

The  next  step  in  this  investigation  was  incorporating  uncertainty  quantification  into  the 
design  optimization.  Gradient-based  reliability-based  design  optimization  incorporating  evidence 
theory  was  utilized.  It  was  discovered  that  in  reliability-based  design  problems,  gradient 
information  is  extremely  computationally  expensive.  The  computational  cost  arises  from  finding 
the  gradients  where  the  number  of  simulations  needed  in  the  reliability  analysis  is  the  number  of 
variables  plus  one.  An  adaptive  particle  swarm  optimization  program  was  developed  to  reduce 
this  computational  time. 

Finally,  three  separate  analyses  were  completed  on  the  RLV  wing  incorporating 
uncertainties.  The  first  was  an  aleatory  uncertainty  quantification  incorporating  structural 
uncertainties  and  an  aeroelastic  flutter  limit  state.  This  analysis  discovered  that  the  FORM  and 
SORM  methods  produced  results  with  nearly  the  same  reliabilities.  This  indicates  that  the  limit 
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state  is  nearly  linear  when  only  structural  parameters  are  used  in  an  aleatory  investigation.  The 
second  analysis  conducted  was  an  evidence  theory  uncertainty  quantification  analysis,  where  the 
uncertainties  were  found  in  atmospheric  conditions  and  composite  materials.  For  this  study  six 
uncertain  variables  were  used,  three  of  which  came  from  structural  uncertainties  and  three  from 
atmospheric  uncertainties.  Flutter  dynamic  pressure  was  used  as  the  limit  state.  This 
investigation  concluded  that  when  epistemic  uncertainties  were  considered  a  wider  reliability 
bound  was  obtained.  A  wider  bound  implied  the  need  for  additional  infonnation  about  uncertain 
parameters.  Although  a  high  reliability  was  determined,  it  was  still  necessary  to  perfonn  a  risk 
minimization  optimization.  The  third  analysis  conducted  was  a  gradient-based  reliability  design 
optimization  of  the  RLV  wing  including  aeroelastic  uncertainties.  The  uncertain  variables  were 
three  sections  of  skin  thicknesses  as  well  as  air  density.  The  design  variables  in  the  optimization 
were  the  three  skin  sections  thicknesses,  where  the  objective  function  was  to  minimize  weight. 
The  optimization  of  the  RLV  wing  demonstrated  that  a  reduction  in  weight  as  well  as  an  increase 
in  reliability  could  be  obtained  through  optimal  modification  of  thickness  distributions. 
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